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1  Motivation  and  Objectives 

One  of  the  most  promising  techniques  for  the  prediction  of  turbulent  flows  is  that  of  Large  Eddy 
Simulation  (LES),  in  which  an  under-resolved  representation  of  the  turbulence  is  simulated  numer¬ 
ically  by  modeling  the  effects  of  the  unresolved  small-scales  on  the  simulation.  Such  simulations 
have  been  applied  in  several  flows  with  reasonable  success.  However,  there  are  several  outstand¬ 
ing  problems  that  need  to  be  addressed  before  LES  can  fulfill  its  promise  as  a  tool  for  turbulence 
prediction  in  engineering  flows.  The  most  serious  problems  limiting  the  usefulness  of  LES  is  the 
representation  of  turbulence  near  walls  and  other  strong  inhomogeneities  and  the  dependence  of 
models  on  the  filter  and/or  numerical  discretization. 

The  optimal  LES  formulation1,2  provides  a  rigorous  framework  in  which  to  address  these  issues 
and  to  develop  and  analyze  LES  models  and  simulations.  Optimal  LES  modeling  has  been  found 
to  produce  accurate  LES  simulations  when  based  on  reliable  statistical  information,  so  the  pri¬ 
mary  thrust  of  the  current  research  is  to  reduce  or  eliminate  the  need  for  empirical  statistical  input 
through  theory  and  modeling  of  turbulence  statistics.  When  small-scale  isotropy  is  a  valid  assump¬ 
tion,  the  Kolmogorov  theory  and  isotropy  are  can  provide  much  information.  However,  when  in¬ 
homogeneity  and  anisotropy  are  strong,  or  the  Reynolds  number  is  not  too  large,  more  information 
will  be  required,  and  models  for  this  are  being  developed,  particularly  for  near-wall  turbulence. 

Theoretical  models  for  the  turbulence  multi-point  correlations  allow  optimal  LES  models  to  be 
implemented  relatively  simply  in  production  CFD  codes,  and  preliminary  implementations  in 
FDL3DI  at  AFRL  have  been  pursued. 
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2  Background  &  Approach 


The  starting  point  for  the  development  of  LES  is  the  definition  of  a  spatial  filter  :,  which  can  be 
applied  to  the  Navier-stokes  equations  to  obtain  an  equation  for  the  filtered  velocity  up. 

diii  duiUj  dp  1  d  2Ui 

~dt  =  ~  dxj  ~  ~dxi  +  Re  dxjdx,  +  () 

Where  Mt  is  the  sub-grid  model  (force)  term,  which  includes  the  divergence  of  the  sub-grid  stress 
as  well  as  terms  that  arise  when  the  filter  does  not  commute  with  differentiation.  The  problem  in 
LES  of  course  is  to  model  M*.  A  very  important  result  of  our  research 1  is  that  an  LES  w  will  match 
the  one-time  statistics  of  filtered  turbulence  u  if  and  only  if  the  model  nrii(w)  of  M%  is  given  by 

m,i(w)  =  ( Mi(u)\u  =  w)  (2) 

This  model  also  minimizes  the  difference  between  M*  and  m  ,  (in  the  mean-square  sense),  and  so 
this  model  has  all  the  properties  that  one  could  ask  of  a  sub-grid  model.  We  therefore  call  it  the 
ideal  sub-grid  model. 

Unfortunately,  the  conditional  average  in  (2)  cannot  practically  be  determined,  since  the  conditions 
are  that  the  entire  filtered  velocity  field  match  the  entire  LES  field.  However,  it  can  be  estimated 
using  stochastic  estimation3  which  is  a  well-established  technique  for  estimating  conditional  av¬ 
erages.  The  result  is  a  class  of  estimation  based  LES  models  as  first  proposed  by  Adrian4.  To 
perform  stochastic  estimation  of  the  model  term  Mt,  one  must  provide  correlations  of  the  LES 
variables  w  with  the  model  term,  as  well  as  correlations  of  LES  variables  with  themselves.  Previ¬ 
ously,  before  the  current  effort,  this  data  had  been  obtained  from  direct  numerical  simulations,  but 
this  is  clearly  unacceptable  for  a  predictive  LES  model.  In  this  research,  we  develop  the  theoretical 
models  of  the  required  correlations  to  be  used  in  an  optimal  LES  model. 

As  will  be  shown  in  the  following  section,  it  will  usually  be  sufficient  to  estimate  the  subgrid  model 
term  M  as  linear  in  the  LES  variables,  thus  the  correlations  which  must  be  determined  are: 

(«j(x,)«fc(xw))  and  (Mi(x)uk(x")).  (3) 

We  do  not,  however,  generally  have  theory  for  these  two-point  correlations  of  filtered  and  model 
quantities.  But  we  do  have  theory  for  multi-point  velocity  correlations,  and  as  will  be  shown,  the 
LES  correlations  (3)  can  be  obtained  by  applying  the  filtering  operator  to  the  following  two  and 
three-point  correlations: 


M(r,  s) 

=  (u(r)u(s)) 

(4) 

B(r,  s) 

=  (u(r)u(s)u(s)) 

(5) 

T(r,s,t) 

=  (u(r)u(s)u(t)> 

(6) 

Thus,  models  for  these  correlation  tensors  are  needed.  When  the  small  scales  are  isotropic,  Kol¬ 
mogorov  inertial  range  theory  applies,  and  this  will  be  sufficient  to  determine  all  but  T,  and  an 
extended  inertial-range  model  to  determine  T  has  also  been  developed.  Furthermore,  correlations 
of  filtered  velocities  with  themselves  can  be  computed  directly  from  an  LES,  thus  it  is  possible  to 
determine  some  of  the  necessary  correlations  dynamically  in  a  running  LES.  In  this  case,  a  model 
is  only  needed  for  the  two-point  third  order  correlations  B. 
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When  the  assumption  of  small-scale  isotropy  and  homogeneity  are  not  valid,  the  models  for  the 
required  correlations  must  be  more  complex.  The  most  common  example  of  this  situation  is  near¬ 
wall  turbulence.  In  the  log-layer  of  wall-bounded  turbulence,  however,  self- similarity  of  the  corre¬ 
lations5  can  be  used  to  represent  the  inhomogeneity.  To  model  anisotropy,  we  extend  the  descrip¬ 
tion  of  single-point  anisotropy  in  terms  of  structure  tensors  6  to  the  two-point  correlation. 

With  models  of  the  multi-point  correlations  in  hand,  and  making  use  of  dynamic  procedures,  prac¬ 
tical  optimal  models  that  do  not  require  DNS  data  can  be  formulated.  And  these  have  been  tested 
in  non- wall-bounded  flows. 


3  Supported  Research 


To  pursue  the  objectives  defined  above,  a  number  of  research  activities  were  pursued  under  the 
current  grant.  These  are  described  briefly  below  and  in  more  detail  in  the  following  subsections, 
and  the  referenced  publications. 


1.  Development  and  Testing  of  Theoretical/Dynamic  OLES  Models:  To  produce  OLES 
models  for  turbulence  that  is  isotropic  in  the  small  scales,  theoretical  expressions  for  the  cor¬ 
relations  M  and  B,  are  readily  available  from  Kolmogorov  inertial  range  theory.  However,  the 
three-point  third-order  correlation  is  more  difficult.  To  avoid  needing  a  priori  information 
on  the  three-point  correlation,  a  dynamic  approach  in  which  correlation  of  LES  variables 
with  themselves  is  determined  from  the  LES  as  it  runs,  was  formulated.  It  was  implemented 
along  with  a  the  theoretical  model  for  B,  and  tested  in  isotropic  turbulence  and  a  free-shear 
layer.  An  initial  implementation  in  a  (compressible)  production  CFD  code  (FDL3DI)  was 
also  pursued. 

2.  Development  of  a  Model  for  the  Three-Point  Correlation:  No  previous  model  for  the 
three-point  third-order  correlation  had  previously  been  developed  because  of  the  complexity 
of  this  object.  A  model  of  this  tensor  was  developed  however,  by  recasting  a  classical  result7 
for  the  most  general  tensor  form,  and  selecting  for  the  scalar  functions  appearing  in  the 
form,  the  simplest  expressions  that  can  possibly  be  consistent  with  the  known  inertial-range 
scaling.  This  produced  a  tensor  representation  with  4  adjustable  constants,  which  was  fit  to 
DNS  data. 

3.  Modeling  the  Two-Point  Correlation  in  Wall-Bounded  Turbulence:  To  formulate  near¬ 
wall  OLES  models,  the  two-point  second  and  third  order  correlations  valid  in  the  log-layer 
are  needed.  We  concentrate  first  on  the  second-order  correlation,  but  the  third-order  correla¬ 
tion  can  be  determined  in  terms  of  the  second-order.  There  are  three  complications  that  have 
been  addressed:  inhomogeneity,  anisotropy  and  finite  Reynolds  number.  Inhomogeneity 
and  finite  Reynolds  number  effects  are  treated  together,  because  in  near-wall  turbulence  the 
effective  Reynolds  number  of  the  turbulence  varies  in  the  wall-normal  direction.  A  general 
anisotropy  model  for  the  two-point  correlation  was  developed  in  terms  of  the  structure-tensor 
formulation  devised  by  Kassinos,  Reynolds  &  Rogers6. 

4.  Testing  and  Evaluation  of  the  Filtered  Wall  Formulation  for  the  Wall-Layer:  In  previous 
work,  the  filtered  wall  LES  formulation  was  devised,  which  avoids  many  of  the  problems  of 
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very  near-wall  turbulence,  essentially  by  filtering  them  out.  This  modeling  approach  was 
further  tested  and  modified  to  characterize  the  performance  of  the  model.  Further  it  is  being 
reformulated  for  finite  volumes,  instead  of  the  spectral  methods  originally  used. 
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4  Development  and  Testing  of 

Theoretical/Dynamic  OLES  Models 


It  was  shown  previously8  that  finite- volume  “filtering”  in  which  discrete  volume  averages  con¬ 
stitute  the  LES  variables,  yields  a  good  approximation  to  the  ideal  LES  of  isotropic  turbulence. 
However,  this  assessment  relied  on  the  use  of  DNS  data  to  compute  the  necessary  correlations. 
The  dependence  of  the  correlations  on  DNS  data  does  not  allow  the  use  of  OFVLES  models  in 
other  flows.  Alternative  methods  for  determining  the  correlations,  which  do  not  rely  on  detailed 
empirical  data,  must  be  developed.  In  this  chapter,  theoretical  and  dynamical  approaches  for  de¬ 
termining  the  correlations  in  isotropic  turbulence  OFVLES  models  are  described. 


4.1  Theoretical  Modeling  of  Correlations 

Theoretical  modeling  of  the  correlations  in  OFVLES  will  be  restricted  to  the  ones  required  for 
the  convective  fluxes  (Mconv)  for  two  reasons.  The  theory  used  in  modeling  the  correlations,  Kol¬ 
mogorov’s  scaling  laws  for  velocity  structure  functions9"13,  is  strictly  valid  only  in  the  limit  of 
an  infinite  Reynolds  number.  Also,  the  viscous  terms’  contributions  to  the  evolution  are  small  far 
from  solid  boundaries.  Close  to  solid  boundaries,  the  current  theoretical  formulations  are  not  valid, 
efforts  to  address  this  problem  are  described  in  section  6. 

In  this  section,  the  convective  fluxes  are  defined  as 

Mij  =  /  UjUknk  dx.  (7) 

J, 9% 

where  the  extra  i  and  j  subscripts  are  introduced  to  represent  all  convective  fluxes  through  a  single 
expression.  With  this  notation,  the  estimation  equation  may  be  rewritten  as 


(' rriijEk )  =  ( MijEk ) . 


(8) 


Substitution  of  the  expressions  for  the  model  and  events  in  the  estimation  equation,  Eq.  (8),  yields 
a  linear  system  of  equations  for  the  model  coefficients: 

{M,,}  =  A, j  +  V  L%k  K>  +  E  (<»? ) .  (9) 

a  a,/3 

(WmMij)  =  AH  (wl)  +  L«k  (W%wl)  +  Qtfkl  ( WkWl  Wl)i  0°) 

o;  a, (3 

(WmWnMij)  =  Aij  (w^wi)  +  Ltjk  Qifkl  ( WkWl  WmWn) •  C1 1) 

a  a,/3 
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To  close  these  equations,  the  following  correlations  of  the  filtered  velocity  must  be  computed: 


ig(a)  =  (My),  (12) 

4m(S>7)  =  (!3) 

Ikm(an)  =  (WkWm)  ,  (14) 

43!m(«.A7)  =  ,  (15) 

(S,  7,  <*)  =  (WmWnMij  >  ,  (16) 

4,>.  P, 7,  <5)  =  (w>W)  >  (17) 


where  s  is  the  surface  over  which  the  flux  is  computed.  Note  that  (w)  =  0,  due  to  homogeneity 
and  Galilean  invariance  of  the  filtered  equation. 


The  integrals  representing  the  effect  of  the  finite-volume  filter  on  the  unfiltered  velocities  can  be 
commuted  with  the  averaging  operator,  yielding  functions  that  are  integrals  of  correlations  of  the 
unfiltered  velocities: 


w 

4fe(s> a) 

4(«>0) 

4 fc(a>  Pi  7) 

4  ki(s,a,P) 


Iijki(a>P>'r>s) 


Wj(x)Mj(x))  dx, 

rq(x)'Uj(x)'Ufc(x/))  dxdx', 

^(x^^x'))  dxdx', 

'uj(x)MJ  (x/)Mfc(x"))  dxdx'dx", 
ui(x.)uj(x)uk('x.l)ui(x"))  dxdx'dx", 

Ui  (x.)uj  (x')-Ufc (x")wn (xw) )  dxdx'dx"  dx" 
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I  ^  J  f3  J  a. 


(3  J  a  J  s 


(18) 

(19) 

(20) 
(21) 
(22) 
(23) 


J'y  J  [3  J  a 

Therefore,  to  close  the  system  of  estimation  equations  for  the  quadratic  OFVLES  model  coeffi¬ 
cients  for  the  convective  fluxes,  the  following  correlations  of  the  unfiltered  velocities  are  required: 


Rij( r)  =  (ui(x)uj(x  +  r)) ,  (24) 

Tijk( 0,  r)  =  («i(x)«j(x)«fc(x  +  r)) ,  (25) 

Tijk{rl,  r2)  =  (ui(x)uj(x  +  r1)wfc(x  +  r2))  ,  (26) 

Fijki{ 0,  r2,  r3)  =  (iii(x)«i(x)«fc(x  +  r2)«z(x  +  r3))  ,  (27) 

Fijkii*1,  r2,  r3)  =  («i(x)«j(x  +  r1)wfc(x  +  r2)tq(x  +  r3))  .  (28) 


Hence,  expressions  for  the  two-point,  second-  and  third-order,  the  three-point,  third-order  and  the 
three-  and  four-point,  fourth-order  correlations  of  the  unfiltered  velocities  are  needed. 


To  model  these  velocity  correlations,  three  assumptions  are  necessary: 


1.  The  turbulence  is  isotropic,  at  least  at  the  length  scales  encountered  in  the  correlations. 
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2.  The  separations  (r,  r1,  r2  and  r3)  are  within  the  Kolmogorov  inertial  range. 

3.  The  Reynolds  number  of  the  turbulence  is  sufficiently  large  that  it  may  be  considered  infinite. 


These  assumptions  are  consistent  with  the  ones  used  by  Kolmogorov  to  deduce  expressions  for 
second-  and  third-order  longitudinal  structure  functions,  which  are  used  to  derive  the  two-point, 
second-  and  third-order  velocity  correlations: 


Rij(  r) 
Tijk{  0,r) 


u25ij  +  ^-e2/3\r\  4/3  (r*'r 3  -  4M%)  , 


(29) 

(30) 


For  the  fourth-order  velocity  correlations,  the  controversial  quasi- normal  approximation  (QNA)  is 
used.  In  this  approximation,  the  fourth-order  cumulants  are  assumed  to  vanish  as  with  normally 
distributed  variables.  This  assumption  allows  the  expression  of  the  fourth-order  correlations  in 
terms  of  products  of  second-order  correlations: 

Fijki( r\  r2,  r3)  =  %  (r1)  Rdj  (r3  -  r2)  +  RVJ  (r2)  Rdj  (r3  -  r1) 

+  Rij  (r3)  Rij  (r2  -  r1)  .  (31) 


While  experimental 14  and  numerical15  evidence  has  shown  that  this  approximation  is  quite  accu¬ 
rate  in  isotropic  turbulence,  use  of  the  QNA  in  two-point  closure  models  has  led  to  negative  energy 
spectra — a  review  of  the  issues  may  be  found  in  Lesieur16.  However,  in  this  case,  the  anomalous 
behavior  of  the  spectra  was  due  to  unrestrained  growth  of  the  third-order  moments  during  the  sim¬ 
ulation  and,  in  the  application  to  OFVLES,  third-order  correlations  are  not  determined  in  this  way. 
It  must  be  emphasized  that  the  moment  equation  simulations  using  QNA  were  quite  accurate  for 
small  time  intervals  and  displayed  unphysical  behavior  only  after  long  simulation  times 16 . 

The  only  remaining  correlation  needed  for  closure  of  the  estimation  equations  is  the  three-point, 
third-order  velocity  correlation.  The  general  expression  for  this  correlation  is  more  easily  ex¬ 
pressed  in  wave-space 14.  The  six-dimensional  Fourier  transform  of  T  and  its  inverse  are  defined 
by 


fijk(  k\k2) 

Tijk(  r4,r2) 


(2tt) 


Tijk{ r1,  r2)e-i(kl  r  +k  r  Wdr 


tijk{ r\r2)ei(k  r  +k  r  MkMk2, 


k1  7  k2 


(32) 

(33) 


and  the  general  expression  for  T  is 


Tijk( k1,  k2)  =  k1  -  k2)'Pj„(k1)'Pfcp(k2)(<I>(fi;1,  k2,  k3)^^  -  $(n3,  n2,  Kl)k2n8mp 

+  <F( K \  K3,  K2)kp5nm  +  'T( K 1,  K2,  K3)klmk2nkl ),  (34) 


7 


where  ft1  =  j kl\,  k2 

defined  by 


| A;2 1  and  ft3  =  \kl  +  k2\.  The  divergence-free  projection  operator  V  is 


’Pij(k)  —  Sij 


k|2  ’ 


(35) 


Using  the  QNA  to  close  the  evolution  equation  for  the  third-order  moments,  Proudman  &  Reid 17 
derived  the  following  equations  for  <f>  and  T : 

J^(«\  «2,  «3)  =  /(«3)  |7V)  -  /(«)]  -  v  (fa1)2  +  (O2  +  («3)2)  $(k\  k2,  ft3)  (36) 
ft2,  ft3)  =  -V  ((ft1)2  +  (ft2)2  +  (ft3)2)  ^(ft1,  ft2,  ft3),  (37) 

where  /(ft)  =  U(k)/47tk2  and  E(k)  is  the  three-dimensional  energy  spectrum. 


For  stationary  turbulence,  the  function  T  vanishes.  For  <f>,  the  eddy-damping  approximation 
(EDQNM)  is  introduced  to  avoid  the  realizability  issues  encountered  when  using  the  QNA  for 
closure.  The  damping  term  actually  reproduces  the  damping  effects  of  the  fourth-order  cumu- 
lants  that  were  discarded  in  the  QNA.  Using  the  EDQNM  approximation  and  assuming  an  infinite 
Reynolds  number  (so  that  the  viscous  term  may  be  dropped)  yield  the  following  equation  for  <f>: 


^(fi;1,  k2,  k3) 


f(n3)  [/(k2)  -  /(«); 


H  (/c1,  ft2,  ft3) 


(38) 


Finding  the  inverse  Fourier  transform  of  <f>  is  equivalent  to  finding  the  real-space  representation  of 
the  three-point,  third-order  correlation  T.  Unfortunately,  the  analytical  inverse  Fourier  transform 
of  <f>  is  complicated  due  to  the  singular  behavior  of  the  integrand  at  the  origin.  There  has  been  no 
solution  to  this  problem  to  date. 

However,  Eq.  (38)  can  still  produce  valuable  information.  The  damping  term  is  usually  of  the  form 

^(ft1,  ft2,  ft3)  ~  e1/,377(ft1,  ft2,  ft3),  (39) 

where  T  is  a  function  only  of  geometric  quantities,  not  of  any  flow-dependent  quantities.  Substi¬ 
tuting  the  expressions  for  the  spectra 


E(k)  =  Ce2/3k~5/3 


(40) 


and  the  damping  term  into  Eq.  (38)  yields  a  function  for  T  of  the  form 

T^ft1,  ft2,  ft3)  =  e7fjifc(ft1,  ft2,  ft3).  (41) 


The  scaling  equation  for  the  Fourier  coefficients  of  the  three-point,  third-order  velocity  correlation, 
Eq.  (41),  may  also  be  found  in  a  more  straightforward  way  by  invoking  Kolmogorov’s  scaling  ar¬ 
guments  for  structure  functions  in  the  inertial  range910.  Kolmogorov  postulates  that,  in  the  inertial 
range  of  turbulence,  the  velocity  structure  function  may  depend  only  on  the  dissipation  and  sepa¬ 
ration,  since  it  is  far  removed  in  scale  from  the  production  and  dissipation  ranges  (where  viscosity 
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and  the  length  scales  of  the  forcing  would  be  of  importance).  Therefore,  by  using  dimensional 
analysis,  clearly  the  only  possible  functional  form  for  the  three-point,  third-order  velocity  correla¬ 
tion  is  linear  in  dissipation,  yielding  Eq.  (41). 

Note  that  the  three-point,  third-order  velocity  correlation  is  a  linear  function  of  dissipation,  e,  mul¬ 
tiplied  by  a  term  is  a  function  only  of  geometrical  or  spatial  properties  of  the  grid/filter,  Tt.  There¬ 
fore,  the  function  for  the  correlation  can  be  divided  into  two  parts,  one  which  depends  only  on  the 
flow,  and  is  a  known  function  of  dissipation,  and  one  that  depends  only  on  geometric  properties  of 
the  numerical  grid/filter.  The  geometry-dependent  function  is  unknown  and  can  be  approximated 
only  through  a  difficult  inverse  Fourier  transform.  However,  for  a  given  grid  configuration,  a  single 
DNS  may  be  used  to  supply  the  correlation  data  for  any  flow  situation. 


Even  in  the  simple  case  of  a  cubic  grid,  such  as  the  one  used  in  the  isotropic  turbulence  simulations 
in  this  dissertation,  it  is  not  possible  to  derive  closed  forms  for  the  integrated  correlations,  with  the 
exception  of  the  two-point,  third-order  correlation,  I1.  Function  / 1  can  easily  be  computed  for  the 
cubic  grid  due  to  the  linear  nature  of  the  integrand: 


e 
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{SikXj  +  Sjkxf )  ^  , 


(42) 


where  x"  is  the  center-point  of  the  cubic  volume  a  and  A  is  the  grid  spacing. 


The  integrals  dependent  on  the  two-point,  second-order  correlations — I2,  J4  and  J5 —  are  singular 
at  the  origin,  r,  while  the  integrand  for  J3  is  not  even  known.  In  this  dissertation,  Gauss  quadrature 
was  used  to  compute  the  integrals  with  known  integrands,  including  I1. 


4.1.1  Theoretical  Scaling  of  the  Estimation  Equations 


The  theoretical  expressions  derived  in  the  previous  section  can  be  used  to  suggest  the  appropriate 
scaling  for  the  integral  functions  I  and  the  coefficients  A,  L,  and  Q 18.  The  scalings  for  func¬ 
tions  and  coefficients  for  the  convective  flux  OFVLES  model,  with  dependencies  on  volumes  and 
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surfaces  removed  to  simplify  notation,  are 


/?.  =  u2A25-I°- 
I1  =  eA6/1- 

l3  l3 

t2  _  2/3  A  20/3  T2 
1ij  t  LA  1  ij 

4  =  eAl04 
4  =  “4A84 
4  =  «4Al24 

Aij  =  u2A25ij 

/  e  y/3  - 


(43) 


Q  ij  ^4  Qiji 


where  the  quantities  with  an  overbar  are  order-one  scaled  quantities. 


Substituting  these  scalings  in  the  estimation  equations,  Eq.  (9)-(l  1),  results  in  the  scaled  estimation 
equations 


4  ~  Aij+  (^u^)  EE ^4’  (44) 

4  =  E 4-4  A  E  E  44’  (45) 

7J  =  (^) 2/3  VS-  +  (§) 4/3  E  V7 + E  E  gvl-  <46> 

Two  length  scales  appear  in  the  estimation  equations:  the  filter  scale,  A,  and  the  large  turbulence 
scale,  u3  / e.  The  ratio  of  the  filter  scale  to  the  large  turbulence  scale  appears  in  Eq.  (44),  multiplying 
the  quadratic  coefficient,  and  in  (46),  multiplying  both  the  constant  and  linear  coefficients.  In  LES, 
the  filter  scale  should  be  in  the  inertial  range,  so  that  the  ratio  of  the  filter  scale  to  the  large- 
turbulence  scale  should  be  small  and,  therefore,  Eqs.  (44)  and  (46)  should  uncouple  from  Eq.  (45). 
In  the  limit  when  this  ratio  tends  to  zero,  A  <<  u3/e,  the  system  of  estimation  equations  becomes 


T 0  —  A.. 

1ij 

(47) 

^-“EVS  +  EEvT 

(48) 

J4  =  V  V  0  P 

1IJ  /  J  /  J  y^WWJ-WWJ 

(49) 

where  Eq.  (47)  is  an  uncoupled  equation  for  the  constant  coefficient  A,  Eq.  (49)  is  an  uncoupled 
equation  for  the  quadratic  coefficient  Q,  and  Eq.  (48)  is  the  equation  for  the  linear  coefficient  L, 
which  is  dependent  on  the  solution  for  Q. 
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4.1.2  Results  for  Theory-Derived  Convective  Flux  Models 


The  comparison  of  theory-  and  DNS-based  OFVLES  models  is  difficult  due  to  the  low  Reynolds 
number  of  the  DNS  simulations  and  the  lack  of  viscous-flux  models  for  infinite-Reynolds-number 
theory -based  OFVLES.  However,  it  is  of  interest  to  compare  the  results  of  the  theory-  and  DNS- 
based  OFVLES  models  even  though  the  theory-based  model  is  used  in  simulations  outside  its 
range  of  applicability.  Closure  of  the  theory-based  models  is  achieved  using  scaled  DNS  data  for 
the  third-order,  three-point  velocity  correlation.  The  values  for  the  velocity  variance,  u2,  and  the 
dissipation,  e,  used  in  computing  the  theory-based  correlations,  were  set  from  the  DNS  data  and 
the  energy  input  for  the  simulation,  respectively. 

A  posteriori  statistics  for  theory-  and  DNS-based  OFVLES,  and  dynamic  Smagorinsky  models  are 
compared  in  Fig.  1.  Note  that  both  the  DNS-based  OFVLES  and  dynamic  Smagorinsky  models 
contain  viscous  terms,  but  the  theory-based  model  does  not.  While  the  theory-based  model  most 
accurately  simulates  the  third-order  structure  function,  especially  at  the  smallest  separations,  the 
energy  spectrum  results  are  not  as  good.  The  large  wavenumber  energy  spectrum  from  the  theory- 
based  model  is  considerably  lower  than  the  filtered  DNS  result  and  that  of  the  other  models.  The 
relatively  low  Reynolds  number  of  the  simulated  flow  may  be  the  reason  for  such  discrepancy, 
given  that  the  theory-based  model  is  strictly  valid  only  in  the  infinite  Reynolds  number  limit. 

A  study  of  the  effects  of  stencil  geometry,  has  been  conducted  for  the  theory-derived  OFVLES 
models,  to  confirm  that  the  conclusions  reached  in  DNS-based  OFVLES 819  are  valid  for  the  theo¬ 
retically  derived  models.  In  addition,  a  grid/filter  refinement  analysis  for  a  high-Reynolds-number 
isotropic-turbulence  case  was  performed,  which  was  not  possible  with  the  DNS-based  OFVLES 
models,  which  are  restricted  to  a  single  Reynolds  number  and  grid  configuration. 


Analysis  of  linear  part  of  the  theory-derived  models  As  in  the  analysis  of  the  DNS -based 
OFVLES  models,  two  quadratic  flux  models  with  identical  quadratic  parts  (with  events  identical 
to  the  S4  DNS-based  model)  but  with  different  sized  linear  stencils  displayed  little  variation  in  both 
a  priori  and  a  posteriori  results.  The  a  priori  variation,  displayed  in  Table  1,  changes  insignifi¬ 
cantly  when  increasing  the  size  of  the  linear  stencil,  consistent  with  the  conclusion  reached  when 
analyzing  the  DNS-based  models.  The  a  posteriori  statistics,  Fig.  2,  also  display  little  variation 
in  the  results  between  the  model  using  the  smaller  1x1x2  stencil  and  the  one  using  the  larger 
1x1x4  stencil.  A  small  difference  in  the  spectra  is  noticed  only  at  the  highest  wavenumbers, 
while  the  structure  functions  for  the  two  models  show  only  minimal  differences  throughout  the 
separation  range. 


With  the  simplification  afforded  by  the  small  stencil  for  the  linear  part  of  the  OFVLES  model  and 
the  decoupling  of  the  quadratic  part,  it  is  possible  to  write  a  simplified  expression  for  the  coefficient 
L  for  the  linear  part: 


Lij 


{ ' MijWi  )  Sa,/3  Qij 


a,  f- 3  — 

W^WjWi 


{wfwi  )  -  {wt  Wi  ) 


(50) 


where  there  is  no  sum  on  repeated  indices.  The  +  and  —  superscripts  denote  velocity  events  in 
volumes  shifted  in  the  positive  and  negative  directions,  respectively,  normal  to  the  flux  face.  The 
numerator  of  Eq.  (50)  represents  the  energy  transfer  due  to  the  linear  part  of  the  model  and  is 
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(b)  Third-order  structure  function,  S3 

Figure  1:  Filtered  DNS  statistics  compared  with  theory-based  OFVLES,  DNS-based  OFVLES  and 
dynamic  Smagorinsky  results. 
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Table  1 :  A  priori  relative  estimation  variation  for  theory -based  models  with  varying  linear  stencil 
geometry  only.  DNS-based  results  repeated  for  comparison. 


Linear  stencil 

DNS-based 

MUW  M ww 

Theory-based 

M-UW  M-WW 

1x1x2 

1x1x4 

0.3165  0.2071 
0.3165  0.2070 

0.3193  0.2179 
0.3193  0.2186 

Table  2:  A  priori  relative  estimation  variation  for  theory-based  models  with  differing  constraints 
for  the  quadratic  part.  DNS-based  results  are  repeated  for  reference. 


DNS-based 

Theory-based 

Quadratic  Stencil 

Flux  Model 

Flux  Model 

Size 

Type 

Muw 

-MWw 

Muw 

-M-ujw 

1x1x2 

Optimal 

0.3165 

— 

0.3193 

— 

Standard 

0.3317 

— 

0.3320 

— 

Optimal 

— 

0.2145 

— 

0.2179 

1x1x4 

2  DOF 

— 

0.2168 

— 

0.2187 

1  DOF 

— 

0.2168 

— 

0.2181 

Standard 

— 

0.2372 

— 

0.2395 

simply  the  theoretical  energy  transfer  of  the  flux  minus  the  energy  transfer  due  to  the  quadratic 
part  of  the  model. 


Analysis  of  quadratic  part  of  the  theory-derived  models  To  corroborate  the  conclusion  arived 
at  previously,  based  on  DNS  data,  that  standard  discretizations  for  the  quadratic  part  of  the  con¬ 
vective  flux  models  are  accurate  representations  of  the  ideal  LES,  with  negligible  difference  in 
accuracy  when  compared  with  optimal  models,  is  repeated  for  the  theory-based  OFVLES  models. 

In  Table  2  the  a  priori  variation  for  the  theory-based  OFVLES  models  are  displayed,  along  with 
DNS -based  variations  for  comparison.  The  trend  for  the  theory -based  models  is  similar  to  the 
DNS-based  ones.  The  a  priori  variation  for  models  with  a  single  DOF  is  only  slightly  larger  than 
that  of  the  optimal  model,  while  the  standard  discretizations  have  variations  that  are  around  10% 
larger.  Notice  that  the  a  priori  variations  of  the  theory -based  models  are  only  slightly  larger  than 
the  DNS-based  ones,  despite  the  fact  that  the  variations  are  computed  based  on  low-Reynolds- 
number  data. 

In  Fig.  3,  a  posteriori  statistics  for  the  theory-based  models  with  varying  quadratic  parts  are  dis¬ 
played.  The  difference  among  the  energy  spectra  of  the  different  models  is  negligible  in  most  of 
the  wavenumber  range.  Only  at  the  highest  wavenumber  is  it  possible  to  distinguish  the  statistics 
of  the  different  models.  The  result  is  similar  for  the  third-order  structure  function,  reinforcing 
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(b)  Third-order  structure  function,  S3 

Figure  2:  Filtered  DNS  statistics  compared  with  theory-based  OFVLES  results  for  varying  linear 
part  stencil  sizes. 
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the  conclusion  reached  previously,  that  is,  that  the  standard  discretization  for  the  quadratic  part  of 
the  convective  flux  model  is  sufficiently  accurate,  being  negligibly  less  accurate  than  the  optimal 
model. 


Grid/Filter  refinement  for  high-Reynolds-number  turbulence  A  high-Reynolds-number  flow 
was  also  simulated  using  theory-based  models.  The  results  obtained  for  this  case  are  compared 
with  a  filtered  infinite-Reynolds-number  theoretical  energy  spectrum  and  a  grid/filter  refinement 
study  was  also  conducted  (Fig.  4). 

No  grid  convergence  is  observed,  since  the  filter  width  is  varying  among  the  different  simula¬ 
tions,  but  the  dependence  of  the  solutions  on  the  filter  width  is  as  expected.  The  three-dimensional 
energy  spectra  of  the  more  refined  models  clearly  becomes  asymptotic  to  the  filtered  theoretical 
energy  spectrum.  The  third-order  structure  functions  converge  to  a  linear  in  r  dependence  in  the 
central  range  of  separation  values,  as  clearly  seen  in  the  1283  simulation,  which  is  representative 
of  inertial-range  energy  transfer.  The  larger  separations  are  affected  by  the  forcing;  the  small  sep¬ 
aration  energy  transfer  is  damped  by  the  filtering.  These  results  demonstrate  that  the  theoretical 
modeling  of  the  velocity  correlations  necessary  for  OFVLES  allows  the  extension  of  optimal  mod¬ 
els  to  high-Reynolds-number  turbulence,  though  a  treatment  for  third-order,  three-point  correlation 
is  needed. 


4.2  Dynamic  Computation  of  Correlations 

In  the  previous  section,  theory-based  and  scaled  DNS  velocity  correlations  were  used  to  extend 
OFVLES  models  to  high-Reynolds-number  turbulence.  Dependence  on  DNS  data  was  restricted 
only  to  the  three-point,  third-order  velocity  correlations  in  a  given  grid  configuration.  To  remove 
the  dependence  on  DNS  data  completely,  an  alternate  way  of  computing  the  three-point,  third-order 
velocity  correlations  is  needed.  It  is  also  important  to  derive  correlations  for  turbulence  that  departs 
from  the  assumptions  made  in  the  preceding  section,  such  as  near-wall  and  low-Reynolds-number 
flows,  to  extend  OFVLES  to  more  complex  flows. 

The  estimate-event  correlations,  (■ mEi ),  appearing  on  the  right-hand  side  of  the  estimation  equa¬ 
tions,  Eqs.  (9)-(l  1),  can  be  rewritten  as  G'?  ( E,  E} ) ,  where  ( E^Ej )  is  a  matrix  of  event-event  corre¬ 
lations,  consisting  only  of  functions  of  the  filtered  velocity,  and  Cj  is  a  vector  of  model  coefficients. 
It  is  possible  to  compute  the  event-event  correlation  matrix  on  the  fly  using  the  velocity  field  from 
an  ongoing  simulation,  thus  eliminating  the  need  of  a  theoretical  model  for  all  terms  on  the  right- 
hand  side  of  Eqs.  (9)— ( 11).  This  dynamic  procedure  for  computing  correlations  is  extremely  ad¬ 
vantageous  in  complex  turbulence,  since  it  limits  the  necessary  theoretical  modeling  to  two-point, 
third-order  and  three-point,  fourth-order  velocity  correlations  for  quadratic  event  estimates  and  to 
only  the  two-point,  third-order  velocity  correlation  when  using  a  standard  discretization  for  the 
quadratic  part  of  the  flux  model. 

A  more  fundamental  reason  for  computing  event-event  correlations  using  a  dynamic  procedure 
rather  than  theoretical  or  DNS-based  expressions  is  the  following.  When  using  theoretical  or  DNS- 
based  expressions  for  the  event-event  correlations,  the  model  coefficients  are  adjusted  a  priori  to 
match  the  estimate-event  correlations  with  the  modeled-term-event  correlations.  On  the  other 
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(b)  Third-order  structure  function,  S3 

Figure  3:  Filtered  DNS  statistics  compared  with  theory -based  OFVLES  results  for  varying  con¬ 
straints  on  the  quadratic  part  of  the  model. 
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(a)  Three-dimensional  energy  spectra,  E{k) 


Figure  4:  Theoretical  three-dimensional  energy  spectrum  and  third-order  structure  function  com 
pared  with  theory -based  OFVLES  results.  Numbers  in  the  legend  represent  grid  points  per  direc 
tion. 
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hand,  the  dynamic  procedure  matches  the  a  posteriori  estimate-event  correlations  with  the  theo¬ 
retical  flux-event  correlations,  as  the  simulation  is  running.  It  can  be  argued  that,  if  done  correctly, 
the  dynamic  procedure  does  not  contain  any  modeling  error  in  the  event-event  correlations,  since 
they  are  a  true  representation  of  the  statistics  in  the  simulation,  not  an  a  priori  estimate. 

On  the  other  hand,  the  dynamic  approach  for  collecting  event-event  correlations  is  limited  by  the 
sample  size  defined  by  the  discretization  and  flow.  Non-homogeneous  grids  and  flows  limit  the 
application  of  space-averaging.  The  small  sample  size  in  these  cases  is  an  issue  in  other  LES 
models,  such  as  the  dynamic  Smagorinsky  model,  and  no  definite  solution  has  been  found  to  date. 


4.2.1  Issues  in  the  Dynamic  Computation  of  Correlations 


Unlike  DNS-  and  theory-based  OFVLES  models,  it  is  not  possible  to  use  all  the  DOF  of  the 
quadratic  part  of  the  flux  model  for  the  optimization  of  the  variance  of  the  model.  After  long 
simulation  times,  the  results  from  the  dynamic  models  exhibit  inaccurate  large-scale  behavior.  It  is 
expected  that  fully  converged  velocity  correlations,  such  as  those  used  in  DNS-  and  theory-based 
simulations,  yield  optimal  models  that  are  compatible  with  the  large-scale  dynamics,  but  poorly 
converged  correlations,  such  as  those  used  during  the  dynamic  procedure,  do  not.  Adding  to  the 
issue,  time-averaging  of  the  three-point,  third-order  correlations  is  exceedingly  difficult  because  of 
long-time  scales  associated  with  the  production-range  (small-wavenumber)  velocity.  These  long¬ 
time  scales  lead  to  slow  convergence  of  the  dynamically  computed  correlations  and  of  the  resulting 
optimal  model.  It  is,  therefore,  important  to  constrain  the  large-scale  behavior  of  the  OFVLES 
dynamic  models  while  the  simulation  is  underway.  To  constrain  the  OFVLES  models  properly,  it 
is  important  to  realize  that  the  large-scale  dynamics  of  the  OFVLES  model  should  be  similar  to 
the  standard  model  derived  from  order-of-accuracy  considerations. 

The  use  of  1  x  1  x2  stencils  for  the  linear  part  is  also  essential  for  the  stability  of  dynamic  OFVLES 
convective-flux  models,  especially  the  Muw  flux  model.  When  using  larger  stencils  for  the  linear 
part,  the  linear  part  of  the  Muw  model  becomes  numerically  (linearly)  unstable.  Note  that,  as 
verified  previously819  and  in  4.1.2,  in  the  cases  of  DNS-  and  theory -based  OFVLES  models,  the 
1x1x2  stencil  was  optimal  for  the  linear  part  of  the  models.  The  addition  of  extra  dependencies 
on  further  removed  velocity  events  did  not  increase  the  accuracy  of  the  models,  a  priori  or  a 
posteriori ,  even  when  accurately  represented.  The  addition  of  these  events  to  the  dynamic  models 
not  only  failed  to  increase  the  accuracy  of  the  model,  as  was  verified  for  the  Mww  flux,  but  caused 
the  model  to  become  unstable,  as  noted  for  the  Muw  flux. 

The  more  pressing  issue  in  the  OFVLES  of  bulk  flow  turbulence  is  the  accurate  computation  of  the 
energy  dissipation  e.  In  the  case  where  the  quadratic  part  of  the  model  is  constrained  to  the  stan¬ 
dard  discretization,  the  only  quantity  that  must  be  computed  using  Kolmogorov  theory  is 
which  is  a  linear  function  of  e,  as  given  by  Eq.  (30).  For  isotropic  turbulence,  the  average  dissi¬ 
pation  can  be  computed  by  using  the  two-point  correlation,  Eq.  (29),  but  for  the  case  of  a  generic 
turbulence  the  problem  has  not  been  solved. 
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4.2.2  Results  for  Dynamic  OFVLES 


The  results  obtained  with  the  dynamic  OFVLES  models  differ  from  the  DNS-  and  theory-based 
results  in  a  fundamental  manner.  In  both  the  DNS-  and  theory-based  cases,  the  values  of  the  dissi¬ 
pation  and  velocity  variance  from  the  DNS  simulations  were  used  to  compute  the  OFVLES  models 
before  the  simulation  was  initiated  (implicitly  in  the  DNS -based  case  and  explicitly  in  the  theory- 
based  case).  For  dynamic  OFVLES,  these  values  must  be  computed  on-the-fly  since  it  is  necessary 
that  the  correlation  values  computed  dynamically  and  theoretically  be  consistent.  Therefore,  at 
every  time  step,  the  dissipation  is  approximated  using  the  two-point  correlation,  Eq.  (29),  and  the 
velocity  variance  is  approximated  by  the  resolved  velocity  variance.  These  approximations  are 
then  inserted  into  the  expressions  for  the  left-hand  side  of  Eqs.  (9)-(10). 

The  statistics  obtained  in  simulations  using  the  dynamic  OFVLES  model  for  the  S2  and  S4  stencil 
are  shown  in  Figs.  5  and  6.  The  number  of  DOF  available  for  the  optimization  of  the  quadratic  part 
was  varied  for  each  of  the  stencils.  Note  that  one  DOF  is  for  the  Muw  flux,  while  the  remaining 
DOF  are  for  the  Mww  flux. 

It  is  clear,  for  both  the  S2  and  the  S4  stencils,  that  unconstrained  dynamic  OFVLES  models  are 
not  accurate  for  the  energy  spectra  throughout  the  wavenumber  range.  The  only  case  where  the 
unconstrained  model  produced  an  accurate  energy  spectrum,  the  S4  stencil  with  3  DOF,  displays 
the  worst  results  for  the  third-order  structure  function.  The  best  results  for  both  the  energy  spectrum 
and  third-order  structure  functions  were  obtained  with  either  a  single  DOF  for  optimization  of  the 
quadratic  part  or  the  standard  discretization,  as  for  the  theory-based  OFVLES  models. 

In  Figs.  7  and  8,  the  dynamic  OFVLES  statistics  are  compared  with  those  of  DNS-  and  theory- 
based  OFVLES  models  and  the  dynamic  Smagorinsky  model,  where  all  models  use  the  S4  stencil. 
In  Fig.  7,  all  models  use  the  standard  descretization  for  the  quadratic  part  of  the  convective  flux, 
while  in  Fig.  8  the  OFVLES  models  use  an  optimal  quadratic  part  with  one  DOF.  Also,  it  is 
important  to  remember  that  both  the  theory-based  and  the  dynamic  OFVLES  models  do  not  contain 
viscous  terms,  while  the  DNS-based  and  dynamic  Smagorinsky  simulation  do.  The  lack  of  viscous 
terms  for  the  OFVLES  models  could  account  for  the  small  differences  between  the  models  at  the 
large  wavenumbers  and  smallest  scales. 

It  is  difficult  to  determine  whether  the  use  of  the  standard  discretization  or  the  optimal  model 
with  a  single  DOF  is  more  accurate.  While  the  third-order  structure  function  results  for  standard 
discretization  seem  to  be  more  accurate  than  those  for  the  models  with  a  one  DOF  optimal  quadratic 
part,  the  opposite  is  true  for  the  energy  spectrum.  However,  these  results  show  that  the  accuracy  of 
the  dynamic  OFVLES  model  is  comparable  with  that  of  other  approaches,  with  the  only  noticeable 
discrepancy  being  the  results  for  the  third-order  structure  function  at  the  smallest  separation. 


4.3  Application:  Simulation  of  a  Temporally-Evolving  Shear  Layer 

In  the  previous  section,  the  optimal  finite-volume  LES  model  was  applied  only  to  isotropic  tur¬ 
bulence.  To  evaluate  OFVLES  in  a  more  complex  flow,  the  simulation  of  a  temporally-evolving 
shear  layer  is  pursued  here.  The  shear  layer  is  appropriate  as  the  next  test  case  given  that  the  flow 
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(b)  Third-order  structure  function,  S3 

Figure  5:  Filtered  DNS  statistics  compared  with  dynamic  OFVLES  results.  The  quadratic  parts  of 
the  OFVLES  models  are  constrained  to  differing  degrees.  S2  stencil. 
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(b)  Third-order  structure  function,  S3 

Figure  6:  Filtered  DNS  statistics  compared  with  dynamic  OFVLES  results.  The  quadratic  parts  of 
the  OFVLES  model  are  constrained  to  differing  degrees.  S4  stencil. 
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(b)  Third-order  structure  function,  S3 

Figure  7:  Filtered  DNS  statistics  compared  with  theory -based,  DNS-based  and  dynamic  OFVLES 
and  dynamic  Smagorinsky  results.  All  models  are  constrained  to  using  Taylor-series  derived  ap¬ 
proximations  for  the  quadratic  part  of  the  convective  flux.  S4  stencil. 
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(b)  Third-order  structure  function,  S3 

Figure  8:  Filtered  DNS  statistics  compared  with  theory-based,  DNS-based  and  dynamic  OFVLES 
and  dynamic  Smagorinsky  results.  OFVLES  models  are  constrained  to  using  1  DOF  for  the 
quadratic  part  of  the  convective  flux.  Dynamic  Smagorinsky  model  uses  the  standard  discretization 
for  the  quadratic  part.  S4  stencil. 
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physics  are  more  complex  without  the  necessity  of  modeling  near-wall  turbulence.  One  can  expect 
that  the  results  obtained  using  the  dynamic  OFVLES  model  in  the  simulation  of  a  shear  layer  is 
representative  of  the  performance  of  OFVLES  models  away  from  the  wall. 


4.3.1  Issues  of  Dynamic  OFVLES  for  the  Shear  Layer 


Compared  with  the  isotropic  turbulence  case,  an  important  difference  in  the  formulation  of  the 
dynamic  OFVLES  model  in  the  shear  layer  is  the  existence  of  a  mean  velocity.  For  shear-layer 
turbulence,  the  estimation  equations  for  the  constant  and  linear  term  must  be  rewritten  as 

(Ma)  =  Aij  +  L±  (wf)  +  Lf  (wf)  +  <%> ,  (51) 

(Mijw;  >  =  Aij  (wf )  +  Lf  (wf wf  )  +  Lf  (wf wf  >  +  (fl’l3wf >  ,  (52) 

(MHwf )  =  Aij  (wf  >  +  Lf  (wf wf  >  +  Lf  (wf  wf  >  +  (%jwf)  ,  (53) 


where 

*ij  =  Y,Qvw?w?  (54) 

a,/3 

is  used  to  simplify  the  notation. 


For  the  current  simulation,  the  linear  coefficients  are  constrained  to  be  of  equal  value  but  opposite 
sign,  Lf  =  —Lf,  leading  to  a  purely  dissipative  model  with  no  dispersion.  With  this  simplifica¬ 
tion,  Eqs.  (5 1)— (53)  become 

(Mij)  =  +  Lij  (wf  -w~)  +  (^ij) ,  (55) 

(Mij  (wf  -w~))  =  A^  (wf  ~w~)  +  Lij  (  (wf  -w~)2^  +  (%j  (wf  -wf)).  (56) 

The  OFVLES  flux  models  on  faces  in  the  homogeneous  directions,  that  is,  with  face  normals 
pointing  in  the  x  and  z  directions,  may  be  simplified  by  noting  that  (wf)  =  (wf),  decoupling 
Eq.  (55)  from  Eq.  (56).  Also  note  that  the  value  of  Atl  is  irrelevant  for  the  momentum  computation, 
since  it  is  constant  for  opposing  faces  in  the  homogeneous  directions  of  the  shear  flow.  Therefore, 
for  the  x  and  z  directions,  the  equation  for  the  linear  term  of  the  dynamic  OFVLES  model  is 

(M^  (wf  -wf))  =  L^  ( (wf  -wf)2)-  (Vij  (wf  -wf)),  (57) 


where  the  solution  for  Lt]  is 


Lij  — 


(M^  (wf  -  Wj))~  (Vjj  (wf  -Wi)) 
(wf-wfY 


(58) 


Note  that,  since  in  the  homogeneous  directions  the  mean  velocities  (wt)  and  (wi  )  are  identical, 
the  equation  for  the  OFVLES  model  is  identical  to  the  one  used  in  homogeneous  turbulence. 


The  estimation  of  the  OFVLES  flux  models  for  faces  in  the  y  direction  is  considerably  more  com¬ 
plicated  due  to  two  factors.  First,  the  value  of  the  constant  Al3  does  contribute  to  the  filtered 
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evolution  equation,  since  flux  models  on  opposing  faces  in  the  y  direction  will  have  different  val¬ 
ues  for  the  constant  coefficient.  Second,  the  equations  for  the  linear  and  constant  parts  of  the  flux 
model  are  coupled,  since  (wt+)  ^  (wj  ). 

The  solution  for  the  OFVLES  model  fluxes  in  the  y  direction  is  given  by 

(Mij)  -  (vl//'i)  -  Lij  < K+  -  «>r)  >  ,  (59) 

((My  -  (Mij))  K+  -  ^r))  -  ((gjj  -  (%))  K+  -  wr))  (60) 

( K+  -  wr)‘ 2)  -  (wf  -  w~ )2 


While  the  dynamically  computed  terms  in  Eqs.  (58)  and  (59)-(60)  are  easily  solved  for  during 
the  simulation,  the  same  is  not  true  for  the  correlations  depending  on  M,  which  are  computed  via 
the  Kolmogorov  relations.  The  correlations  developed  using  Kolmogorov’s  inertial-range  theory 
do  not  account  for  the  existence  of  a  mean  velocity,  being  correct  only  for  the  correlations  of 
the  velocity  fluctuations.  However,  as  is  shown  in  Zandonade19,  if  small-scale  isotropy  of  the 
velocity  correlations  is  a  valid  approximation  in  shear  flow,  the  theoretical  expressions  for  these 
correlations,  developed  from  Kolmogorov’s  inertial  range  theory,  may  still  be  used. 


To  compute  the  third-order,  three-point  velocity  correlations  using  the  theoretical  expressions  de¬ 
veloped  above,  it  is  necessary  to  estimate  the  dissipation  at  the  face  where  the  flux  is  computed. 
For  the  shear-layer  simulation,  the  dissipation  is  approximated  using  the  scaling  equation 

<  W  -  wi)  W  -  wJ )  >  =  ct/'^  (61) 

where  the  scaling  coefficient  Ce  is  computed  using  the  Kolmogorov  relation  for  the  second-order 
structure  function 


=^L 
ij  3V2 


-  Y-|-4/3 


V  Jv 


(x+ —  x  (x+ —  x  )  .  —  4|x+ —  x  \25ij  j  dx  dx+ 


-  Y”|-4/3 


JV+ Jv~ 

Therefore,  the  dissipation  is  computed  using  the  relation 


x+  —  x  ) .  (x+  —  x  )  .  —  4|x+  —  x  \25ij  j  dx  dx+  j  .  (62) 


e  = 


((<  ^  )  (wt~wj)) 


3/2 


c: 


(63) 


The  dissipation  is  computed  for  every  face  using  the  appropriate  velocity  for  a  given  grid  and  then 
averaged  in  each  xz  plane. 


4.3.2  Details  of  the  Simulation 


The  geometry  for  the  simulation  of  the  shear  layer  is  displayed  in  Fig.  9,  where  the  x  and  z  direc¬ 
tions  are  periodic.  The  initial  condition  is  imposed  by  setting  the  plane-averaged  mean  velocity  in 
the  x  direction  to 


<«>**  (y)  =  AU  erf  (f )  > 


(64) 
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Figure  9:  Geometry  of  the  shear  layer  simulation. 


where  6  is  a  parameter  setting  the  initial  thickness  of  the  shear  layer,  AU  is  the  velocity  difference 
across  the  shear  layer  and  the  subscript  on  the  averaging  operator  defines  the  plane  over  which  the 
averaging  takes  place.  The  velocity  fluctuations  are  computed  with  a  random  number  generator 
using  a  normal  distribution  and  variance  given  by 


2/  n  2  (  y2\ 

a  (: y)  =  q  exp  I  — —  I  , 

where  q 2  is  the  velocity  variance  at  the  origin  of  the  y  coordinate  system. 


(65) 


The  filtered  evolution  equation  is  advanced  in  time  using  the  standard  fourth-order  Runge-Kutta 
scheme.  The  dynamic  OFVLES  model  is  used  for  the  convective  fluxes,  with  a  standard  second- 
order  divergence-form  finite-volume  model  for  the  quadratic  terms.  The  linear  part  of  the  OFVLES 
model  depends  on  the  two  velocity  events  neighboring  the  flux  face,  a  stencil  shown  to  be  optimal 
in  previous  chapters.  The  correlations  needed  for  the  OFVLES  model  are  averaged  in  each  xz 
plane,  and  are  computed  at  every  step  of  the  Runge-Kutta  scheme.  The  viscous  fluxes  are  not 
computed  in  this  simulation,  so  that  the  Reynolds  number  of  the  simulated  flow  is  infinite. 


Values  of  pressure  are  determined  using  the  classical  fractional  step  method  of  Chorin20.  The 
solution  of  the  filtered  evolution  equation  is  a  three-step  process.  First,  the  momentum  equation 
is  solved  without  the  pressure  term,  yielding  a  tentative  solution  for  the  velocity  at  a  time  t  +  5t. 
Second,  a  Poisson  equation  for  pressure  is  solved  that  guarantees  that  the  velocity  at  the  end  of  the 
iteration  will  be  solenoidal.  Finally,  the  pressure  gradient  is  subtracted  from  the  tentative  velocity 
to  yield  the  final  solenoidal  velocity  field: 


w*  =  At  A  jMij  + 


d2w]1 


At 


d2p  dw* 


dxiXi  dxi 


w?+1  =  w* 


=  H, 

aA 

OX i 


Re  dxjXj  ' 


(66) 

(67) 

(68) 


where  H  is  the  divergence  of  the  tentative  velocity  solution.  Note  that  Eq.  (67)  is  simply  the 
divergence  of  Eq.  (68),  where  dw™+1/dxi  is  set  to  zero. 
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The  periodicity  of  the  x  and  z  directions  allows  the  use  of  the  Fourier  transform  to  reduce  the  cost 
of  solving  the  Poisson  equation.  Fourier  transforming  Eq.  (67)  in  the  x  and  z  directions  yields: 


At 


(69) 


The  second  derivative  of  pressure  in  the  y  direction  can  be  computed  by  a  second-order  finite- 
difference  method: 


P(x 


A, 


A2 

y 


A2 

^ y 


P(x)  + 


p(x~ 

A2 

y 


H 
A ? 


(70) 


where  Ay  is  the  grid  separation  in  the  y  direction.  Equation  (70)  is  a  tridiagonal  system  of  equations 
that  is  solved  using  the  tridiagonal  matrix  algorithm  (TDMA)  for  each  kx  and  kz  wavenumber, 
except  for  the  (0,  0)  mode,  which  is  set  to  zero.  The  boundary  conditions  in  the  y  direction  are 
Neumann,  where  the  derivative  of  pressure  in  the  y  direction  is  set  to  0. 


Also  note  that  the  kx  and  kz  wavenumbers  appearing  in  Eq.  (70)  are  substituted  by  effective 
wavenumbers  for  the  second-order-accurate  first-derivative  finite-difference  operator.  With  this 
substitution,  the  current  scheme  is  equivalent  to  a  standard  second-order  finite-difference  dis¬ 
cretization  of  Eq.  (67). 


Periodic  boundary  conditions  for  the  velocity  are  imposed  at  the  x  and  z  boundaries.  For  the 
boundaries  in  the  y  direction,  irrotational  flow  is  imposed.  The  vanishing  vorticity  at  the  boundary 
is  imposed  by  using  the  y  derivative  of  the  continuity  condition: 

d2u  d2v  d2w 
dxdy  dy2  dzdy 


components  of  vorticity  are  zero: 


yielding  a  Poisson  equation  for  v: 


can 

be  expressed 

as  derivatives  of  v  given  that  the  x 

and  z 

dw 

dv 

(72) 

^ X 

=  0  -> 

dy 

dz’ 

du 

dv 

(73) 

Ulz 

=  0  -»• 

dy 

dx' 

V2v 

=  0. 

(74) 

The  periodicity  of  the  domain  allows  the  Fourier  transformation  of  Eq.  (74)  in  the  x  and  z  direc¬ 
tions: 

(i  2  .  /  2\  -  /  /  /  \  .  9v  (y,  kx,  kz) 

~  {kx  +  kz)v(y,kx,kz)  + - —2 - =  0.  (75) 

Substituting  a  second-order  finite-difference  approximation  for  the  derivative  in  Eq.  (75)  yields  a 
difference  equation  for  the  discrete  velocity: 


Vj+i  ( kx ,  kg)  -  (2  +  k2)  vj  (kx,  kz)  +  vj- 1  (kx,  kz) 


A2 

y 


=  0, 


(76) 
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where  k2  =  k2  +  k2  and  the  subscripts  containing  j  are  the  y  indices  of  the  grid  points. 

The  difference  equation  (76)  may  be  solved  by  noting  that  the  solution  must  be  of  the  form: 

Vj  =  rj.  (77) 


Substituting  Eq.  (77)  into  Eq.  (76)  yields: 


Vj  (kx,  kz) 


±  VkPAy 


vo  (kx,  kz) , 


(78) 


where  vq  (kx.  kz)  is  the  Fourier  coefficient  for  an  interior  xz  plane.  To  guarantee  that  v  vanishes  as 

j  ->  cx),  Vj  (. kx ,  kz)  is  given  by 


Vj  ( kx .  kz ) 


Vk?  Ay 


vo  (■ kx ,  /cz) . 


(79) 


With  v  at  the  boundary  set  by  Eq.  (79),  the  corresponding  boundary  conditions  for  u  and  w  may  be 
found  by  using  Eqs.  (72)  and  (73)  with  first-order  finite-difference  approximations  for  the  deriva¬ 
tives. 


4.3.3  Results  for  the  Shear  Layer  Simulation 


The  temporally-evolving  shear  layer  displays  self-similar  growth  while  the  growth  rate  is  linear.  In 
the  linear-growth-rate  regime,  quantities  should  be  scaled  by  the  velocity  difference  A U  and  the 
momentum  thickness 

2^ 


5m  = 


(«>a 


A  U 


d  y. 


(80) 


The  initial  momentum  thickness  is  denoted  by  5q.  The  appropriately  scaled  time  variable  is 

tAU 


T  = 


(81) 


The  results  in  this  section  were  obtained  using  the  parameters  in  Table  3.  In  that  same  table,  the 
grid  parameters  for  the  OFVLES  are  compared  with  those  used  in  the  DNS  by  Rogers  &  Moser21 . 
This  DNS  will  be  used  as  the  benchmark  when  comparing  the  results  from  the  OFVLES. 

The  growth  of  the  shear  layer  during  the  simulation  is  displayed  in  Fig.  10,  while  the  instantaneous 
growth  rate  is  shown  in  Fig.  11.  The  growth  rate  is  approximately  linear  in  the  range  r  =  200  — 
500.  During  this  time  period  the  growth  of  the  shear  layer  will  be  considered  self-similar  and  the 
statistics  will  be  collected.  Note  that  the  growth  rate  is  slightly  higher  than  the  averaged  value  for 
the  self-similar  portion  of  the  DNS.  The  offset  between  the  start  of  self-similar  growth  in  the  DNS 
and  in  the  LES  is  due  to  the  use  of  different  initial  conditions  in  each  of  the  simulations.  In  the 
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Table  3:  Grid  parameters  for  OFVLES  and  DNS  shear  layer  simulations. 


Parameter 

X 

OFVLES 

y  z 

X 

DNS 

y 

z 

Grid  points  (or  modes) 
Size  of  domain 

128 

12550 

128  32 

8050  31.25Jo 

512 

125(50 

210 

oo 

192 

31.2550 

Figure  10:  Growth  of  shear  layer:  dynamic  OFVLES  compared  with  DNS. 
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Figure  11:  Growth  rate  of  shear  layer:  instantaneous  value  for  dynamic  OFVLES  compared  with 
averaged  value  for  self-similar  portion  of  DNS. 


DNS,  turbulence  profiles  from  a  boundary-layer  DNS  were  used  to  start  the  simulation,  while  in 
the  LES  a  synthetic  field  (described  in  Section  4.3.2)  was  used  as  the  initial  field. 

The  large-scale  structure  of  the  simulated  shear  layer  at  time  r  =  300  may  be  visualized  in  Figs.  12 
and  13,  which  plot  the  instantaneous  vorticity  magnitude  and  velocity  vectors  in  an  xy  plane.  The 
dynamic  OFVLES  results  display  large-scale  features  similar  to  the  ones  depicted  in  the  spanwise 
vorticity  contour  plots  by  Rogers  &  Moser21  for  velocity  fields  in  the  linear  growth  regime.  The 
accurate  description  of  the  large-scale  structures  is  expected,  given  that  the  standard  discretization 
scheme  used  as  the  quadratic  part  of  the  OFVLES  model  is  accurate  for  the  large  scales  and  the 
linear  part  of  the  model  is  most  important  for  scales  close  to  the  filter  scale. 

The  profiles  for  the  mean  velocity  and  Reynolds  stresses,  averaged  during  the  self-similar  portion 
of  the  simulation,  are  shown  in  Figs.  14-18.  The  results  obtained  with  the  OFVLES  model  are 
comparable  to  those  obtained  with  the  DNS,  except  for  the  (vv)  Reynolds  stress,  which  was  higher 
in  the  OFVLES  results.  This  discrepancy  may  be  due  to  the  approximation  of  the  mean  flux  used 
to  compute  the  flux  models  for  the  faces  with  normals  pointing  in  the  y  direction.  These  results 
are  evidence  of  the  robustness  of  the  OFVLES  model  in  the  case  of  bulk-flow  turbulence,  since 
accurate  results  were  obtained  even  though  approximate  expressions  were  used  for  computing  the 
dissipation  and  the  mean  flux. 
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Figure  12:  Shear  layer  visualization:  vorticity  magnitude  in  an  xy  plane. 
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Figure  14:  Mean  velocity  profile  for  shear  layer:  dynamic  OFVLES  compared  with  DNS. 


Figure  15:  Mean  ( uu )  Reynolds  stress  for  shear  layer:  dynamic  OFVLES  compared  with  DNS. 
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Figure  16:  Mean  (vv)  Reynolds  stress  for  shear  layer:  dynamic  OFVLES  compared  with  DNS. 


Figure  17:  Mean  (ww)  Reynolds  stress  for  shear  layer:  dynamic  OFVLES  compared  with  DNS. 
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Figure  18:  Mean  (uv)  Reynolds  stress  for  shear  layer:  dynamic  OFVLES  compared  with  DNS. 


4.4  Conclusions 

The  validity  of  the  stochastic  estimate  as  an  approximation  to  the  ideal  LES  was  previously  tested 8 
by  using  DNS  data  for  the  necessary  correlations.  The  OFVLES  models  developed  from  the  DNS 
data  were  shown  to  be  more  accurate  and  less  expensive  than  current  LES  models  such  as  the 
dynamic  Smagorinsky  model,  even  if  much  more  limited  in  their  applicability  due  to  the  reliance 
on  the  DNS-based  velocity  correlations.  The  important  properties  of  OFVLES  models  for  isotropic 
turbulence  were  uncovered  in  this  phase  of  the  research.  It  was  determined  that  the  optimal  stencil 
of  velocity  events  in  the  linear  part  of  the  stencil  was  that  of  the  two  events  closest  to  the  face 
where  the  flux  is  modeled.  Larger  stencils  led  to  models  which  had  almost  identical  a  priori  and 
a  posteriori  behavior  as  the  smallest  stencil,  at  a  much  greater  cost.  Another  important  feature  of 
OFVLES  models  is  the  fact  that  the  accuracy  of  the  quadratic  part  of  the  model  computed  through 
the  optimal  procedure  is  similar  to  that  of  a  standard  numerical  discretization  for  the  convective 
flux.  Application  of  these  two  constraints  on  OFVLES  models  lead  to  inexpensive,  yet  accurate 
LES  models  for  isotropic  flow. 

The  second  stage  of  the  development  of  OFVLES  models  was  concerned  with  removing  the  de¬ 
pendence  of  the  models  on  DNS  data.  To  this  end,  it  was  necessary  to  develop  alternate  methods 
to  compute  the  velocity  correlations  appearing  in  the  estimation  equations.  Theoretical  models  for 
the  velocity  correlations  are  derived  from  Kolmogorov’s  inertial-range  expressions.  A  purely  the¬ 
oretical  model  has  not  been  developed  at  this  time,  due  to  one  velocity  correlation,  the  three -point, 
third-order  correlation,  remaining  unclosed.  Also,  a  dynamic  approach  for  computing  resolved 
velocity  correlations  is  developed  to  substitute  DNS  data,  including  the  three-point,  third-order  ve¬ 
locity  correlation.  Results  from  both  the  theory  and  dynamic  approaches  for  determining  velocity 
correlations  are  quite  encouraging,  being  only  slightly  less  accurate  than  the  DNS-derived  models. 
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Note  that  the  theory  used  to  determine  the  velocity  correlations  is  valid  only  in  the  limit  of  an 
infinite  Reynolds  number  and  the  DNS  simulations  used  as  a  benchmark  were  at  a  low  Reynolds 
numbers.  The  OFVLES  models  based  on  dynamic  and  theory-based  correlations  do  not  depend  on 
any  DNS  data,  and  are  the  first  optimal  models  to  display  this  characteristic. 

The  application  of  the  OFVLES  models  to  a  more  complex  flow,  the  temporally-growing  shear 
layer,  was  also  successful.  The  results  obtained  for  mean  velocity  and  Reynolds  stresses  in  the 
LES  are  comparable  to  the  DNS  results  of  Rogers  &  Moser21,  while  using  a  low-order-accurate 
numerical  method  and  a  mesh  sixteen  times  coarser  in  the  x-z  plane  than  the  DNS.  These  results 
are  important  because  they  are  an  indication  of  the  performance  of  OFVLES  models  in  bulk-flow 
turbulence. 

Although  OFVLES  has  been  successfully  applied  in  the  simulation  of  isotropic  and  shear-flow 
turbulence,  several  issues  remain  to  be  resolved  before  these  models  can  meet  widespread  use: 


•  Dissipation:  To  evaluate  the  flux-velocity  correlations  that  are  derived  from  theory,  it  is 
necessary  to  compute  the  averaged  dissipation  at  the  face  where  the  flux  is  modeled.  In 
isotropic  turbulence  many  different  approaches  exist  for  computing  the  dissipation,  based 
mainly  on  Kolmogorov’s  expressions  for  the  second-  and  third-order  structure  functions;  the 
second-order  function  was  used  in  this  work.  On  the  other  hand,  in  more  complex  flows  no 
standard  procedure  for  computing  dissipation  has  been  developed.  For  these  flows,  an  accu¬ 
rate  method  for  computing  the  dissipation  is  necessary  in  order  for  the  OFVLES  methods  to 
remain  useful.  The  optimal  model  simply  ensures  the  correct  distribution  of  the  dissipation 
among  the  different  fluxes  and,  therefore,  if  the  dissipation  is  not  computed  accurately,  the 
optimal  model  itself  will  be  in  error. 

•  Averaging  for  the  dynamic  procedure:  In  the  flows  analyzed  in  this  dissertation,  the  ex¬ 
istence  of  homogeneous  directions  in  the  flow  eased  the  task  of  averaging  the  dynamically 
gathered  correlations.  In  more  complex  grid  geometries — note  that  even  a  simple  flow  ge¬ 
ometry  with  a  complex  grid  arrangement  would  suffer  from  this  issue — it  is  not  clear  how 
to  proceed  with  the  averaging  of  the  velocity  correlations.  While  spatial  averaging  may  be 
used  in  the  case  of  isotropic  turbulence  and  other  flows  that  have  the  appropriate  spatial 
symmetries  (such  as  plane-by-plane  averaging  for  plane  channel  flow),  for  complex  geome¬ 
tries  some  kind  of  time  averaging  must  be  employed.  The  time  average  used  for  computing 
the  correlations  will  affect  the  evolution  of  the  LES,  and  finding  the  correct  procedure  for 
averaging  is  not  trivial.  The  issue  is  similar  to  the  averaging  of  the  dynamic  Smagorinsky 
coefficient,  so  it  may  be  possible  to  use  many  of  the  same  tools  in  the  case  of  OFVLES. 

•  Three-point,  third-order  velocity  correlation:  While  the  dynamic  approach  would  en¬ 
counter  problems  in  flows  with  no  homogeneous  directions  due  to  the  necessity  of  averaging 
the  correlations,  a  purely  theoretical  model  would  not  suffer  from  this  issue.  In  fact,  the  only 
further  work  involved  when  applying  the  theoretical  model  to  a  completely  unstructured  grid 
is  extra  bookkeeping.  A  model  of  the  three-point  third-order  correlation,  which  would  enable 
such  a  theoretical  model,  is  discussed  in  section  5. 

•  Models  for  near-wall  turbulence:  While  all  LES  models  lack  appropriate  wall-modeling, 
the  problem  in  OFVLES  is  focused  on  developing  appropriate  velocity  correlations  near  the 
wall.  Development  of  models  for  near- wall  LES  are  described  in  sections  6  and  7 
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•  Filtering/discretization  in  inhomogeneous  grids:  Inhomogeneous  grids  are  the  norm  in 
complex  or  wall-bounded  geometries.  In  these  grids,  the  varying  filter  scale  yields  an  inho¬ 
mogeneous  representation  of  the  filtered  variables;  that  is,  the  representation  of  the  variables 
is  strongly  dependent  on  the  position  in  the  grid.  An  appropriate  deconvolution  method  is  a 
necessity  for  filtered  simulations,  especially  for  inhomogeneous  grids,  in  order  to  place  the 
results  in  different  simulations  and  even  in  different  regions  of  a  flow  geometry  on  a  common 
ground. 

•  Application  to  other  discretization  schemes:  While  the  work  described  here  has  focused 
on  the  finite-volume  discretization,  the  methods  presented  can  readily  be  applied  to  finite- 
difference  and  finite-element  methods.  The  only  difference  in  the  modeling  procedure  is  a 
redefinition  of  the  filter  and,  therefore,  of  the  quadratures  used  to  compute  the  velocity  cor¬ 
relations.  A  preliminary  study  with  finite-difference  methods  has  been  initiated,  but  a  more 
in-depth  analysis  is  necessary.  No  reason  exists  to  believe  that  any  of  these  discretization 
methods  is  better  or  worse  suited  for  LES  and,  in  order  to  be  widely  accepted,  the  OFVLES 
model  must  be  adapted  to  the  different  discretization  procedures. 
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5  Development  of  a  Model  for  the  Three-Point  Correlation 


Multi-point  velocity  correlations  are  central  to  the  statistical  description  of  homogeneous  isotropic 
turbulence2223.  The  two-point  second-order  velocity  correlation,  and  its  Fourier  transform,  the 
velocity  spectrum  tensor,  are  the  most  commonly  considered  correlations,  and  a  variety  of  models 
for  their  evolution  have  been  developed  (e.g  EDQNM 16). 

In  the  evolution  equation  for  the  two-point  velocity  correlation  in  homogeneous  turbulence,  the 
two  point  third  order  correlation  appears  because  of  the  non-linear  terms  in  the  Navier-Stokes 
equations.  In  isotropic  turbulence,  the  two-point  correlation  equation  reduces  to  the  Karman- 
Howarth  equation.24  In  combination  with  the  Kolmogorov  inertial  range  scaling  assumptions,  this 
leads  to  the  well-known  Kolmogorov  4/5  law25  for  the  third  order  longitudinal  structure  function 
in  the  inertial  range. 

Isotropy  of  the  turbulence  also  implies  that  the  two-point  third-order  correlation  can  be  uniquely 
determined  in  terms  of  the  third-order  structure  funcion,  just  as  the  two-point  second-order  cor¬ 
relation  can  be  writen  in  terms  of  the  second-order  structure  function.  The  two-point  third  order 
correlation  describes  the  transfer  of  energy  from  large-scales  to  small,  and  as  such  is  of  criti¬ 
cal  importance  to  the  theory  of  the  two-point  statistics  of  turbulence.  It  is  also  this  relationship 
to  energy  transfer  that  allows  the  third-order  structure  function  to  be  determined  under  the  Kol¬ 
mogorov  assumptions,  without  an  adjustable  constant. 

Of  course,  the  two-point  third-order  correlation  is  a  restricted  case  of  the  three-point  correlation. 
But  there  is  no  direct  association  of  the  more  general  three-point  quantity  with  the  two-point  cor¬ 
relation  equation  or  with  energy  transfer.  There  has  thus  been  little  motivation  to  study  the  three- 
point  correlation,  and  because  of  the  complexity  of  this  quantity,  there  has  been  virtually  no  work 
to  characterize  it. 

In  Large  Eddy  Simulation  (LES),  however,  there  is  a  motivation  to  characterize  the  three-point 
third  order  correlation.  To  see  why  this  is  true,  note  that  the  evolution  equation  for  the  two-point 
correlation  of  an  LES  field  includes  a  term  involving  the  two-point  third-order  correlation  of  the 
LES  (i.e.  filtered)  velocity  field,  that  arises  from  the  quadratic  terms  in  the  LES  equations.  This 
correlation  of  filtered  velocities  represents  energy  transfer  among  scales  in  the  LES  due  to  the 
nonlinear  term  and  is  thus  of  importance  to  analyzing  the  dynamics  of  the  LES  equations.  It  is 
also  one  of  the  statistical  inputs  to  the  optimal  LES  modeling  approach  1,8  26 .  The  two-point  third 
order  correlation  of  filtered  velocity  can  be  determined  by  applying  the  filter  (three  times)  to  the 
three-point  third-order  correlation  (see  section  5.1.1),  which  is  the  reason  for  our  interest  in  the 
three-point  correlation. 

The  notable  exception  to  the  lack  of  work  on  the  three-point  third-order  correlation  is  the  the  paper 
by  Proudman  &  Reid7,  in  which  the  most  general  isotropic  form  of  the  Fourier  transform  of  the 
three-point  correlation  is  derived.  In  this  paper,  we  start  with  some  mathematical  background  (sec¬ 
tion  5.1).  Then  from  Proudman  &  Reid’s  result  for  the  Fourier  transform  of  the  correlation,  the 
equivalent  form  for  the  physical- space  correlation  is  determined  and  we  find  the  simplest  expres¬ 
sion  consistent  with  the  Kolmogorov  4/5  law  (section  5.2).  Conclusions  and  implications  are 
discussed  in  section  5.3. 
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5.1  Background 


The  Kolmogorov  inertial  range  theory  for  high  Reynolds  number  turbulence  yields  expressions 
for  the  longitudinal  structure  functions,  under  the  well  known  similarity  assumptions  that  in  the 
inertial  range,  the  statistical  properties  of  turbulence  depend  on  the  seperation  scale  and  the  rate  of 
dissipation22-27.  The  longitudinal  structure  functions  are: 

SP(r)  =  <(u||(x)  -  «||(x  +  r))p)  =  Cp(er)p/a  (82) 

where  r  =  |r|  is  the  magnitude  of  the  separation  vector  r,  which  is  assumed  to  be  in  the  inertial 
range,  u\\  is  the  velocity  component  in  the  separation  direction,  e  is  the  average  rate  of  kinetic 
energy  dissipation  (per  unit  mass)  and  Cp  are  the  Kolmogorov  constants,  which  are  generally 
determined  empirically  (e.g.  C-2  ~  2).  The  Kolmogorov  expressions  for  the  structure  functions  are 
found  to  be  quite  accurate  for  p  =  2  and  3,  and  their  accuracy  degrades  as  p.  One  of  the  remarkable 
consequences  of  the  theory  is  that  the  evolution  equation  for  the  two-point  correlation,  in  which 
the  third-order  two-point  correlation  appears,  requires  that  O3  =  —4/5,  the  so-called  Kolmogorov 
4/5  law. 


Isotropy  and  the  continuity  constraints  are  sufficient  to  determine  the  second-  and  third-order  two- 
point  correlation  tensors  from  the  second  and  third  order  structure  functions  respectively.  Using 
the  Kolmogorov  expressions  above,  the  correlation  tensors  are: 


^ijk  (r) 


u2$ij  +  Ye2/3(r)  4/3 (uo  -  4 (r)2Sij) 


where  the  two-point  correlations  are  defined: 

%0)  =  (u(x)ty(x  +  r)) 

%k{  r)  =  (vi(x)vj(x)vk(x  +  r), 


(83) 

(84) 


(85) 

(86) 


and  u2  is  2/3  the  turbulent  kinetic  energy  which  is  also  the  velocity  variance.  The  result  for  the 
second-order  correlation  is  well  known,  but  we  are  not  aware  of  a  previous  reporting  of  the  third- 
order  two-point  correlation  as  shown  here,  though  it  is  implicit  to  the  derivation  of  the  4/5  law,  and 
is  aluded  to  by  Frisch22.  For  completeness,  an  outline  of  the  derivation  is  given  in  Appendix  I. 


The  three-point  third-order  velocity  correlation,  which  is  the  quantity  of  interest  here  is  defined: 

Tijfc(r,  s)  =  (vi(x)vj(x  +  r)vk(x  +  s))  (87) 

The  three  points  form  a  triangle,  and  vectors  r  and  s  determine  it’s  size,  shape,  and  orientation. 
Because  it  is  a  function  of  two  vector  arguments,  T  is  much  more  complex  than  the  two-point 
correlation  S.  However,  we  seak  the  analog  of  (84)  for  the  three-point  correlation.  That  is  a 
“simple”  tensor  form  consistent  with  known  constraints  and  the  Kolmogorov  thoery,  particularly 
the  4/5  law. 
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5.1.1  Relationship  to  LES 


In  the  usual  formualtion  of  Large  Eddy  Simulation  (LES)  the  filtered  velocity  is  defined: 

Uj(x)  =  j  G(x  —  x')vi(x')dx'  (88) 

where  G(x)  is  the  homogeneous  filter  kernel.  The  Navier-Stokes  equations  are  then  filtered  to 
arive  at  an  evolution  equation  for  v: 

d  Vi  dviVj  1  dP  d2v  dr u 

-  _|_  -  — - — J—  U -  — J—  - 

dt  dxj  p  dxi  dxjdxj  dxj 

where 

Tij  =  ViVj  ~  VpJj 

is  the  subgrid  stress.  Our  concern  for  the  moment  is  not  with  the  subgrid  stress  but  with  the 
convection  term. 


(89) 

(90) 


When  deriving  the  evolution  equation  for  the  two-point  correlation  of  the  filtered  velocity  M^-(r)  = 
(vi(x)vj(x  +  r)),  the  convection  term  gives  rise  to  terms  that  can  be  writen  as  derivatives  with 
respect  to  separation  of  the  two-point  third  order  correlation  of  filtered  velocities 

Sijk(r)  =  ( Vi(x)vj(x)vk(x  +  r),  (91) 

in  the  same  way  that  the  unfiltered  two-point  third  order  correlation  arises  in  the  unfiltered  two- 
point  correlation  equation.  As  in  the  unfiltered  equation,  the  energy  transfer  between  scales  of  the 
filtered  velocity  is  mediated  by  this  term,  including  any  transfer  to  the  unresolved  sub-filter  scales. 
Because  §  includes  the  product  of  filtered  velocities  evaluated  at  the  same  point  x,  it  cannot  be 
determined  by  filtering  §.  It  can,  however  be  found  by  filtering  T: 


G'(s)G(s  —  r')G(s  +  r  —  s,)Tjjfc(r',  s')  dsdr'  ds' 


(92) 


which  can  be  derived  easily  by  applying  the  filter  (88)  seperately  to  each  of  the  velocities  in  the 
definition  of  T. 


The  two-point  correlation  equation  also  includes  terms  arising  from  the  sub-filter  stress  term  in  the 
LES  equations,  which  are  derivatives  of  the  two-point  correlation  of  the  sub-filter  stress  with  the 
velocity 

Qijk(r)  =  (rij(x)vk(x  +  r)).  (93) 

The  definition  of  tvj  (90)  means  that  0  can  be  expressed 

Qijfc(r)  =  §ijk( r)  -  Stffc(r)  (94) 

where  Sijk( r)  =  (upfj(x)vk(x  +  r))  can  be  determined  by  filtering  the  two-point  third-order  cor¬ 
relation  of  the  unfiltered  velocity 


r)  =  /  /  G(s)G(s  +  r  -  r/)§iifc(r/)  dsdr1. 


(95) 


40 


In  addition  to  it’s  importance  in  analyzing  the  energy  transfer  in  LES  as  described  above,  §  is  also  a 
necessary  input  for  LES  modeling  using  the  optimal  LES  approach 1S.  In  optimal  LES,  models  for 
the  divergence  of  the  subgrid  stress  are  formulated  to  minimize  mean  square  error  using  stochastic 
estimation,  which  requires  the  correlations  described  above.  Thus  for  the  purposes  of  analysis 
of  LES  and  for  optimal  LES  modeling,  we  are  motivated  to  develop  an  analytic  model  for  the 
three-point  third  order  correlation. 


5.1.2  The  Fourier  Transform  of  T 


Proudman  and  Reid7  determined  a  general  form  for  the  Fourier  transform  of  T  in  both  r  and 
s  (this  is  a  six-dimensional  Fourier  transform).  For  an  incompressible,  homogeneous,  isotropic 
turbulence,  the  most  general  possible  form  for  the  Fourier  transform  $  of  T  is  given  by 

p)  =  Aim(q)Ain(k)Afcp(p)  [5  npkm  0  +  A?ipPn0  1  +  $  mnkp  02  +  kmpnk,pC\  (96) 

where  the  wavevectors  k,  p  and  q  are  interrelated  k  +  p  +  q  =  0,  and  Ami(k)  =  dirn  —  kikm/ k2  is 
the  divergence-free  projector.  The  scalar  functions  0,  0i,  ©2  and  (  depend  only  on  the  magnitudes 
of  the  wavevectors.  For  an  outline  of  the  derivation  of  (96),  see  Appendix  II.  Symmetries  in  the 
tensor  T  imply  symmetries  among  scalar  functions: 

Hk,P,q)  =  —4>{Pi  k,  q)  =  0i  (q,k,p)  =  Mk,q,p)  (97) 


Proudman  &  Reid7  also  analyze  the  dynamic  equation  for  <I>  in  the  context  of  the  quasi-normal 
approximation  to  find  independent  (model)  dynamic  equations  for  0  and  These  equations  imply 
that  for  stationary  turbulence,  Q  is  zero.  We  will  thus  assume  that  (  =  0,  and  with  the  symmetries 
expressed  in  (97),  <I>  is  determined  through  (96)  by  a  single  scalar  function  0  of  k,  p  and  q.  We 
start  with  this  form  in  developing  our  real-space  model  for  T. 


5.2  Inertial-range  model  of  T 

To  construct  a  model  for  the  three-point  third-order  correlation,  a  general  tensor  form  consistent 
with  (96)  is  derived,  and  then  the  scalar  function  appearing  in  the  expression  is  selected  for  con¬ 
sistency  with  the  Kolmogorov  4/5  law. 


5.2.1  A  general  form  for  T  in  real  space 


To  develop  the  analog  of  (96)  in  real  space,  it  will  be  inverse  Fourier  transformed  to  yield  an 
expression  for  T.  However,  to  simplify  the  computations  in  real-space,  it  is  convenient  to  recast 
the  expression  as 

$yfc( k,  p)  =  A,:m(q)Ajn(k)Afcp(p)  [Snpikm4>(k,p,  q)+5mpipn4>(p,  q,  k)+8mnikp(f>(k,  q,p)\  (98) 


41 


where  Ami(k)  =  k25mi— kmki  is  a  modified  divergence  free  operator,  and  <j)(k,p,q)  =  —ubik,  p.  q)/(kpq)'2 
is  a  modifies  scalar  function  that  has  the  same  symmetry  properties  as  q b.  The  advantage  of  this 
form  is  that  the  inverse  Fourier  transform  will  not  give  rise  to  inverse  Laplacian  operators.  An 
inverse  Fourier  transform  of  (98)  yields 

T s)  =  V^V^V^upd^ir,  s,  t)  +  5mpdrJ)(t ,  r,  s)  +  5mndspil)(t ,  s,  r)]  (99) 


which  is  thus  our  general  expression  for  T  in  stationary,  homogeneous,  isotropic  incompressible 
turbulence.  Here,  the  third  seperation  vector  is  t  =  r  —  s,  the  scalar  function  ip(r,  s,t)  is  the 
inverse  Fourier  transform  of  <f>,  and  r,  s  and  t  are  the  magnitudes  of  the  seperation  vectors  r,  s  and 
t  respectively.  The  operators  appearing  in  (99)  are  defined: 


d[ 

% 

dt 


d_ 

dsi  r 

d_ 

9  n  s 
d_ 
dr. 


d_ 

dsi 


(100) 

(101) 

(102) 


vskip  4pt 


V?„  =  SijdgdZ  -  d?d« 


1 3 


(103) 


It  is  straight-forward  to  confirm  that  the  expression  for  T  in  (99)  satisfies  the  relevant  symmetry 
and  continuity  constraints  for  the  third-order  three-point  correlation,  provided  that 


# r,s,t)  = -ip(s,r,t),  (104) 

which  is  the  analogue  of  (97).  The  constraints  on  T  are: 

d\Tijk  =  dy:ijk  =  dkTijk  =  0  (105) 

Tijfc(r,s)  =  Tifcj-(s,r)  (106) 

Tyfc(r,  s)  =  Tifci(-t,  -r)  (107) 

(108) 

The  tensor  form  given  in  (99)  is  clearly  linear  in  i\),  indeed  it  can  be  expressed  as: 

T  ijk  =  Cijki^)  (109) 


where  Cl]k  is  the  tensor-valued  linear  operator  implied  by  (99).  To  complete  the  model  of  the 
three-point  third-order  correlation,  we  need  only  specify  ^(r,  s,  t )  satisfying  (104). 


5.2.2  Scalar  function  ^  in  the  inertial  range 


Our  primary  interest  is  a  model  for  T  that  is  valid  in  the  inertial  range,  analogous  to  the  inertial 
range  expression  for  S  (84).  Kolmogorov’s  4/5  law  constrains  S  to  vary  linearly  with  seperation. 
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Since  T  must  reduce  to  §  when  r,  s  or  t  are  zero,  this  linearity  must  be  reflected  in  T  as  well.  More 
generally,  the  Kolmogorov  similarity  argument22-25  requires  that  in  the  inertial  range 

T(ar,  as)  =  aT(r,  s)  (110) 

The  simplest  way  to  ensure  this  linearity  is  to  choose  to  be  a  polynomial  in  r,  s  and 

t.  Since  each  term  in  (99)  is  a  seventh  derivative  of  ?/),  only  terms  with  total  degree  of  8,  will 
contribute  to  the  linear  scaling  of  T.  This,  along  with  the  symmetry  constraint  on  ^  (104)  suggests 
that  be  constructed  from  terms  of  the  form 

pa,b  =  (rasb  -  rbsa)tc  (111) 

with  a  +  b  +  c  =  8  and  a,  b,  c  >  0  There  are  only  20  expressions  of  this  form,  and  of  these  14 
produce  non-zero  T  when  substituted  for  i/j  in  (99). 

However,  all  of  these  14  non-trivial  T  are  singular  when  r,  s  or  t  are  zero.  For  example,  terms 
such  as:  ryr/rpr / s3  arise,  which  is  clearly  singular  at  s  =  0.  In  addition,  terms  like  5ijSk‘r/s  arise, 
which  is  discontinuous  at  s  =  0.  It  was  found,  however,  that  there  is  a  5-dimensional  null  space  of 
the  singular  and  discontinuous  terms.  There  is  thus  a  5-dimensional  space  of  possible  functions 
that  yield  non-singular,  continuous  T.  The  space  is  spanned  by  the  following  5  functions: 


V4  =  -r^:[-27po,3  -  3p0,5  +  4p2,3  +  18p3,5]  (112) 

5760 

ip2  =  1155840  t~315p°-3  +  4po-7  +  56p2’5  -  14°P3,4  +  1260p3)5]  (113) 

=  n gLnfJ-4^0-1  -  i935po,3  -  40pi  2  +  80pi  3  -  60pi  4 
1257000 

+  16^1,5  +  180^2,3  +  990^3,5]  (114) 

f 4  =  *  n  [— 4P0,1  -  1215p0,3  -  36 Pi  2  +  64p!  3  -  36p,  4  +  4plfi  +  108p2,3  -  20p2,5  +  40p3,4  +  270J+J) 

402 l 20 

V4  =  10684800  ^~60po-1  ~  16065Po,3  ~  504p1[2  +  840pi,3  -  420pM 

+24plj7  +  1260p2,3  +  7560p3l5]  (116) 


Where  pl:J  are  as  defined  in  (1 1 1)  above.  These  functions  have  been  normalized  so  that  each  of  the 
Tn  =  £(-0")  satisfies 

Vjk( 0,  r)  =  ^  (SijTk  -  ^{SikVj  +  6jkri)^  (1 17) 

which  is  just  (84)  with  e  set  to  1.  The  analytic  model  we  seek  for  T  is  thus  given  by: 

5  5 

Tyfc(r,  s)  =  ^2  anl+fc(r,  s)  with  an  =  e  (118) 

72— 1  72—1 


While  the  scaler  basis  functions  wn  are  relatively  simple  to  write  down  (112-116),  the  basis  tensors 
T”  are  not.  Indeed  the  expressions  are  so  complex,  they  will  not  be  writen  out  here.  But  a  program 
is  available  to  evaluate  the  tensor  symbolically  and  numerically  at  http  :  /  /turbulence  .  ices  .  utexas  .  ed 
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ai/e 

0.884 

a2/e 

-2.692 

a3/e 

-6.099 

a  4  j  6 

-5.853 

a5/e 

14.760 

Table  4:  Values  of  the  model  coeficients  found  by  fitting  DNS  data  of  Langford  &  Moser1 


To  display  the  features  of  the  five  basis  tensors  defined  above,  we  examine  the  various  components 
of  the  tensor  for  two  special  arrangements  of  the  seperation  vectors.  First  is  with  the  seperation 
vectors  r  and  s  colinear  (parallel,  designated  by  ||),  which,  without  loss  of  generality,  we  choose  to 
be  in  the  X\  direction  (r  =  re1,  s  =  sex).  In  this  case,  there  are  only  seven  non-zero  components, 
of  which  only  T, , ,  and  t|22  are  independent.  The  other  5  (T212,  t|21,  t|33,  t|13  and  t|31)  are 
related  to  t|22  through  symmetry. 

The  second  seperation  vector  configuration  is  with  r  and  s  orthogonal  (designated  by  _L).  Again, 
without  loss  of  generality  r  =  rei  is  chosen  to  be  in  the  x i -direction,  and  s  =  se2  is  chosen  in 
the  a^-direction.  In  this  configuration,  there  are  14  nonzero  components,  of  which  seven  are  inde¬ 
pendent:  T^n,  T^12,  T|21,  T^22,  T^33,  ¥333,  ¥323.  Each  of  these  is  one  of  a  pair  of  symmetrically 
related  components. 

Since  the  tensor  functions  vary  linearly  with  separation,  the  tensors  can  be  normalized  by  p  = 
max(r,  s,  t ),  which  for  the  special  seperation  configurations  considered  leaves  only  the  dependence 
on  s/r.  The  non-zero,  non-redundant  components  of  T  H/p  and  T ±/p  are  shown  in  figure  19  as  a 
function  of  9  —  arctan(,s/r).  Note  that  the  five  basis  tensors  have  similar  structure,  and  that  two 
of  them  are  quite  similar.  There  has  been  no  effort  to  orthogonalize  the  basis. 


5.2.3  Fitting  to  DNS  data 


To  determine  the  5  coefficients  (ai,  a2, ...,  a5}  in  (118),  a  least-squares  fit  to  data  from  a  Direct 
Numerical  Simulation  (DNS)  of  forced  isotropic  turbulence  at  Re \  =  164 1  is  performed.  Let 
E(r,  s)  =  Tdns  —  TmodeI,  be  the  error  tensor.  Then  the  firing  was  done  to  minimize  the  objective 
function: 

F  =  ^  f  EU(r’  s)EJ,-fc(r,  s)  dr  ds+^  J  E ±k(r,  s)E+fc(r,  s )  dr  ds  (119) 

under  the  constraint  that  )P?)  an  =  e,  where  only  separation  vectors  r  and  s  that  are  parallel  or 
perpendicular  are  considered,  to  reduce  the  data  requirements  to  a  manageable  level,  and  the  inte¬ 
grals  are  taken  over  the  domain  in  r  and  s  for  which  p  is  in  the  approximate  inertial  range  for  the 
DNS  (p  G  [0.15,  0.25]).  This  objective  was  selected  as  a  (crude)  approximation  to  the  integral  over 
all  r  and  s  in  the  inertial  range.  The  coeficients  obtained  from  this  fit  are  given  in  table  4.  The 
coefficient  of  determination  is  R 2  =  0.96,  indicating  that  our  model  describes  the  DNS  data  quite 
well. 

The  ability  of  the  model  to  represent  the  DNS  correlations  is  shown  in  figure  20,  in  which  non- 
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Figure  19:  Basis  functions  for  the  non- zero,  non-redundant  components  of  / p  and  T ±/p  (see 

text  for  definitions)  as  functions  of  6  =  arctan(s/r).  [need  color  key^ 
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Figure  20:  DNS  data  and  analytical  model  for  T 


zero  components  of  T / p  are  ploted  as  a  function  of  6,  for  the  parallel  and  perpendicular  seperation 
vectors,  as  in  figure  19.  The  agreement  between  model  and  correlations  is  very  good.  Further 
indication  of  the  quality  of  the  model  is  given  in  figure  21,  where  contour  plots  show  the  non-zero 
components  of  T  as  functions  of  r  and  s  in  both  the  model  and  the  DNS.  Since  there  is  a  symmetry 
in  each  term  shown,  the  DNS  and  model  are  shown  together  in  each  frame,  with  a  line  of  symmetry 
dividing  them.  The  model  and  DNS  are  very  similar.  But,  there  is  a  minor  discrepancy  for  r  and  s 
near  zero,  which  is  due  to  viscous  effects  not  represented  in  the  model. 


5.3  Discussion  and  Implications 

It  is  remarkable  that  the  simple  considerations  of  isotropy  and  Kolmogorov  scaling  are  sufficient  to 
completely  determine  the  two-point  third-order  correlation,  a  third-ranked  tensor.  The  three-point 
third-order  correlation  is  a  much  more  complicated  object,  so  it  is  equally  remarkable  that  the  same 
scaling  and  isotropy  considerations,  along  with  a  plausable  modeling  ansatz  regarding  functional 
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Figure  21:  DNS  data  and  analytical  model  for  T 
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forms  (1 1 1),  is  sufficient  to  specify  a  model  for  the  three-point  correlation  in  the  inertial  range  with 
just  four  free  constants.  The  model  appears  to  fit  low  Reynolds  number  DNS  data  quite  well.  It 
would  also  be  useful  to  test  the  model  against  higher  Reynolds  number  DNS  data. 

While  the  considerations  leading  to  the  model  are  simple,  the  model  itself  is  algebraically  very 
complex.  A  special-purpose  tensor  algebra  program  was  written  to  perform  the  necessary  manip¬ 
ulations.  Given  the  complexity  of  the  expressions,  it  may  be  that  the  the  ability  to  evaluate  tensor 
components  numerically  will  be  most  useful,  and  a  program  to  perform  such  evaluations  has  been 
developed,  and  is  available  at  http  :  /  / turbulence  .  ices  .  utexas  .  edu. 

The  three -point  third-order  correlation  is  of  particular  interest  in  the  analysis  and  modeling  of  LES, 
because  by  applying  the  LES  filter  to  the  three-point  correlation  one  can  determine  third-order  cor¬ 
relations  of  the  filtered  velocity.  This  is  important  in  analyzing  the  transfer  of  energy  among  scales, 
and  in  formulating  LES  models.  In  particular,  the  representation  of  the  three-point  correlation  will 
allow  optimal  LES  models  of  the  type  evaluated  by  Zandonade  et  al 8  to  be  formulated  without  the 
need  for  DNS  statistical  data. 


Determination  of  §  in  the  inertial  range 

The  derivation  of  (84)  starts  from  the  general  form  of  an  isotropic  two-point  third-order  tensor 
function  of  a  vector  argument 


Sijfe(r)  =  (ui(x)u,-(x)ufc(x  +  r))  =  ar^fk  +  bSjkri  +  cSikrj  +  d6ijrk  (120) 

where  the  scalars  a,  b,  c  and  d  are  functions  of  the  magnitude  of  the  separation  vector  r  =  ||r|| 
only.  Symmetry  in  i  and  j  requires  that  b  =  c.  Further,  the  continuity  constraint  ()Yi]k/drk  =  0 
allows  the  functions  a,  b  and  d  to  be  eliminated  in  terms  of  the  third-order  longitudinal  correlation 
function: 

/(r)  =  (^|(x)vn(x  +  r)),  (121) 

where  vn  is  the  velocity  component  parallel  to  the  separation  vector  r.  The  result  is 

§ufe(r)  =  ^7?  +  ^2  +  6*rj)  ^  /)  -  •  (l22) 


The  third-order  longitudinal  correlation  function  is  directly  related  to  the  third-order  structure  func¬ 
tion,  which  in  the  Kolmogorov  inertial  range  is  Ss(r)  =  —  |er.  The  correlation  /(r)  can  thus  be 
written 


f(r) 


Hr) 

6 


2  er 
~15 


(123) 


Substituting  into  122  then  immediately  yields  (84). 
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The  Most  General  Form  for  $ 


Presented  here  is  a  condensed  version  of  the  derivation  from  Proudman  and  Reid 7 .  The  three-point 
third-order  velocity  correlation  is  T.y-fc(r,  s)  =  (vj(x)vj(x  +  rjiyfx  +  s)).  It’s  Fourier  transform  is 

<hiifc(k,p)  =  i(27r)-6  j  J  T^-fe(r,  s)e_l(k'r+p's)  dr  ds  (124) 

Consistency  with  continuity  requires: 


(k  +  Pi)®ijk{  k,  p)  =  kjQijki  k,  p)  =  pk^ijk(  k,  p)  =  0  (125) 

While  the  most  general  isotropic  third-ranked  tensor  function  of  two  vectors  is: 

0mnp(k,  p)  01  kmknkp  T  02 kmknp.p  T  03 kmPnkp  T  04 

T05 PmPnPp  +  (pfiPmPnkp  T  07 Pm^nPp  T  08 ^mPnPp 

T09 km^np  T  010 kn^mp  T"  011  kpSmn  T  012 Pm&np  T  013 Prifimp  T  014Pp^m(il26) 

where  {0i,  02, 0i4}  are  scalar  functions  of  the  magnitudes  of  the  wavevectors  k  =  ||k||,  p  = 
|| p|| ,  and  q  =  ||q||,  and  k  +  p  +  q  =  0. 

To  enforce  incompressibility,  we  employ  the  divergence-free  projector  Aim(k)  =  <5im  —  kikm/k'2. 
So,  to  satisfy  all  three  incompressibility  conditions  in  equation  (125)  and  isotropy,  we  apply  three 
orthogonal  projectors  to  0mnp,  the  result  is  the  most  general  form  for 

^ijfe(k)  p)  Ajm(q)  Ajn(k)  Afcp(p)0mnp(k,  p)  (127) 

Furthermore,  the  triple  projection  operator  directly  eliminates  all  but  4  components  of  0mnp  shown 
in  equation  (126),  so  effectively  the  above  equation  becomes 

^ijfc(k,  p)  Ajm(q)  Ajn(k)  Apjt(p)  [^03 krripn kp  T  09 km&np  T  011  kpSmn  T  013Pn^»np]  (128) 

Renaming  the  scalar  functions  as  follows:  03  — *  C,  09  — 5 ►  0,  0n  — i ►  02,  013  — >  0i,  we  obtain  (96). 
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6  Modeling  the  Two-Point  Correlation  in 
Wall-Bounded  Turbulence 


For  a  given  velocity  field  u(x),  the  two  point  correlation  tensor  R,j  is  given  by 

Rij(x,  r)  =  (-u'(x)'u' (x  +  r))  (129) 

where  u'  =  u  —  (u)  is  the  fluctuating  velocity  field.  Two  point  correlations  of  velocity  fluctuations 
are  important  in  turbulence  modeling  due  to  the  rich  information  they  contain  about  the  local 
and  global  structure  of  turbulent  flows.  More  importantly  for  the  current  purposes,  multipoint 
correlations  are  used  to  construct  Optimal  LES  (OLES)  models 1,8,28 . 

To  construct  fully  analytical  models  for  LES  we  need  to  be  able  to  reconstruct  the  multipoint 
correlations  of  u'  from  the  statistics  of  the  filtered  fields,  which  are  finite  dimensional  quantities. 
Therefore  it  is  important  that  finite  dimensional  representations  of  the  correlations  of  u'  be  devised 
so  that  a  correspondence  can  be  developed  between  the  LES  statistics  and  the  correlations.  It  is 
this  need  for  building  representations  for  velocity  correlations  in  a  turbulent  flow  with  complex 
geometry  that  motivates  the  current  work. 

Most  of  the  models  for  two  point  correlations  have  been  expressed  in  the  fourier  wave-number 
space  in  terms  of  the  spectrum  tensor  (e.g.  theories  based  on  EDQNM29  and  RDT30).  While 
these  models  capture  many  of  the  fundamental  processes  of  turbulence,  the  global  nature  of  the 
fourier  basis  functions  means  that  the  representations  of  Rij  constructed  from  them  cannot  be  used 
easily  for  flows  in  complex  geometries,  especially  in  wall  bounded  flows.  Therefore  we  choose  to 
construct  our  representation  for  Ri3  in  physical  space,  in  terms  of  tensor  functions  of  x  and  r. 

One  of  the  key  properties  of  R.l3  that  needs  to  be  represented  is  it’s  anisotropy,  both  with  regard 
to  how  different  (x.  r)  is  from  (x.  r)5y/3  for  a  constant  x,r  (componental  anisotropy)  and 
how  Rap{x.,  r)  behaves  over  r/r  for  a  given  x,  a,  (3  and  r  (directional  anisotropy).  In  order  to 
address  these  issues,  Procaccacia31  formulated  a  basis  for  Rrj  in  terms  of  subspaces  which  are 
invariant  to  rigid  rotation  of  the  frame  of  reference  (i.e.  an  S0(3)  decomposition)  where  each  finite 
dimensional  subspace  is  constructed  from  a  particular  spherical  harmonics  mode  1^(0,  6).  It  was 
proposed  that  there  was  a  hierarchy  of  subspaces,  where  the  subspace  with  lower  /  contributes  the 
most  to  Rij( r)  at  intermediate  separations3.  This  approach  could  potentially  yield  finite  dimen¬ 
sional  models  for  RtJ ,  but  it  is  unclear  whether  the  proposed  hierarchy  exists  for  turbulent  flows  in 
complex  geometries.  For  instance,  in  a  wall  bounded  inhomogenous  flow,  the  SO(3)  symmetry  is 
broken  in  the  equations  for  R,n  and  it  becomes  hard  to  argue  that  the  proposed  hierarchy  will  be 
true.  This  points  to  a  need  for,  first,  identifying  the  important  dynamic  and  kinematic  properties  of 
Rij  (which  are  finite  in  number);  second,  ensuring  that  the  finite  dimensional  representation  of  RK) 
is  capable  of  capturing  these  properties  in  any  flow  geometry. 

In  this  work  we  construct  a  representation  for  Rl3  in  terms  of  structure  tensors,  proposed  by  Kassi- 
nos  et  al6  (referred  from  now  on  as  KRR).  These  second  order  tensors  are  single  point  moments 
of  derivatives  of  fluctuating  stream-functions,  which  can  therefore  be  related  to  integrals  of  two 
point  correlations  over  r.  Even  though  these  tensors  are  single  point  moments,  they  contain  infor¬ 
mation  about  the  “structure”  of  turbulence,  e.g.  the  distribution  of  energy  in  different  components 

'The  term  intermediate  separation  will  be  used  throughout  this  work  to  denote  r  in  the  inertial  range  »]<r<  C, 
where  77  is  the  kolmogorov  length  scale  and  C  is  the  scale  deli  ned  by  the  geometry  of  the  fbw 
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of  velocity,  the  dimension  of  turbulence,  etc  (discussed  in  detail  in  section  6.1).  In  the  context  of 
constructing  a  closure  model  for  reynolds  stresses,  it  was  shown  by  KRR  that  it  is  vital  to  represent 
the  pres  sure- strain  correlation  in  terms  of  structure  tensors  in  order  to  obtain  the  correct  evolution 
of  the  reynolds  stress  components  in  the  rapid  distortion  limit.  In  light  of  the  useful  properties  of 
structure  tensors,  we  propose  that  one  should  be  able  to  extract  the  correct  value  of  structure  ten¬ 
sors  from  any  representation  of  Ry .  Since  Ry  is  not  unique  to  a  given  set  of  structure  tensors,  our 
approach  is  to  use  the  theory  of  invariants32  to  formulate  the  most  general  linear  representation  of 
Rij  in  terms  of  the  structure  tensors  (section  6.2.1).  We  then  devise  two  sets  of  constraints;  the  first 
set  ensures  that  Ry  satisfies  the  continuity  equation,  and  the  second  set  ensures  that  Ry  has  the 
necessary  conditions  such  that  the  structure  tensors  calculated  from  the  correlation  are  the  same  as 
the  structure  tensors  used  in  the  representation  (section  6.2.2). 

In  a  similar  approach  Oberlack5  has  represented  Rt  j  (x.  r)  in  terms  of  r  and  a  tensor  length  scale 
related  to  the  Reynolds  stress  tensor  and  the  dissipation  rate  tensor  el3  =  v  The  sig¬ 

nificant  difference  in  our  approach  from  the  latter  is  that  unlike  the  structure  tensors,  Reynolds 
stress  and  dissipation  tensors  do  not  have  information  about  the  directional  anisotropy  of  Ry,  and 
therefore  a  better  representation  of  the  directional  anistropy  can  be  expected  in  our  approach. 

In  order  to  evaluate  the  capabilities  and  shortcomings  of  our  representation,  we  obtain  a  model  for 
Rij  by  fitting  it  to  correlation  obtained  from  DNS  and  comparing  the  two  correlations  (section  6.3). 
This  exercise  points  to  improvements  needed  in  the  model. 


6.1  Background 

6.1.1  Definition  of  structure  tensors 

Structure  tensors  are  single  point  moments  of  derivatives  of  stream-functions.  The  velocity  field 
and  stream-function  are  related  to  each  other  through  the  following  equations: 

u'i  =  eitsip'st  (130) 

tfk,k  =  0  (131) 

Equation  (131)  ensures  that  -0'  is  related  to  the  vorticity  through  Poisson’s  equation: 

Qmn^n,m  ^i,tt  (132) 


In  this  work,  our  discussion  and  use  of  structure  tensors  will  be  limited  to  the  specific  case  of 
homogeneous  turbulence.  The  set  of  structure  tensors  which  are  nonzero  for  homogeneous  tur¬ 
bulence  are  given  by  the  Reynolds  stress  (or  componentality)  By,  dimensionality  Dy,  circulicity 

Fi3  and  strophylosis  Q*]k-  They  are  defined  in  terms  of  the  fluctuating  velocity  stream-function  as 
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follows: 


Bij 

tipqCjt.s  q s,t) 

(133) 

Dij 

=  {^'n,i^'n,j) 

(134) 

i,n*P  j,n) 

(135) 

Q*jk 

(136) 

where  Qijk 

=  -  (Uj^i,k) 

(137) 

and  S  is  an  operator  which  symmetrizes  any  third  order  tensor  as  follows: 


S[Tijk\ 


(138) 


Expanding  eqn  (133),  and  using  eqns  (134),  (135)  and  (131),  the  following  relation  can  be  derived 
among  the  structure  tensors: 


Bij  +  DtJ  +  Fij  —  5ijq 2  (139) 

where  q 2  =  Bkk •  Thus  only  two  of  the  three  second  order  structure  tensors  are  linearly  independent 
for  homogeneous  turbulence,  and  we  choose  B  and  D  for  our  modeling  purposes.  It  was  shown  by 
KRR  that  in  the  presence  of  frame  rotation,  the  strophylosis  Q*  becomes  important,  and  therefore 
we  will  use  it  for  our  model  as  well. 


6.1.2  Exact  relation  between  structure  tensors  and  Rl3 


To  express  structure  tensors  in  terms  of  Ri3,  we  define  a  new  tensor  (r) 
Using  eqn  (132),  we  get 


d2$ijk{r) 
dr  fir  i 


_^/  ^(x  +  r)\ 

drk  >  dr  fin  / 

^(“i(xM(x+r)> 

d2Rjn(  r) 

£ imn 


<«j(x)^U(x  +  r)>- 


(140) 


Assuming  that  the  right  hand  side  of  eqn  (140)  vanishes  at  infinity  (enabling  us  to  use  the  free 
space  green’s  function  of  the  laplacian),  and  noting  that  Ql]k  =  <Tijfc(0),  the  following  is  obtained: 


<S[Qijk\  —  B 


£ imn  r  1  d2Rjn(r )  ' 

An  Jy  r  drkdrm 


(141) 
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Where  V  =  E'!.  The  relation  for  Bij  and  Dtj  can  now  be  given  as  follows 


l 

f  1  d2Rij{r) 

tilmQljm 

4:71  , 

/  -i  n  dr 

Jv  r  orkork 

(142) 

1 

f  1  d2Rkk( r)  , 

^ ilmQmlj 

4n  , 

Jv  r  dridrj 

(143) 

Eqn  (142)  can  be  seen  to  be  the  evaluation  of  inverse  laplacian  of  the  laplacian  of  Rij(r)  at  r  =  0, 
and  thus  represents  the  Reynolds  stress  tensor,  which  simply  gives  the  amount  of  energy  in  different 
velocity  fluctuation  components.  DlJ  is  a  measure  of  how  RtJ  (r)  varies  with  r/r  and  this  will  be 
explained  in  the  following  subsection. 


6.1.3  Information  contained  in  DtJ 


In  order  to  understand  the  link  between  the  anisotropy  of  Rtj(r)  and  Drj ,  let  us  examine  a  hypo¬ 
thetical  (but  physically  legitimate)  correlation  given  as  follows: 


Rij  (r)  =  +  /  (r )  [riTlmrm]  Si:j  ( 1 44) 

where  fir)  has  the  properties  which  ensure  that  the  integral  in  eqn  (142)  does  not  become  singular 
and  that  for  r  <  L,  R,  j  ~  Bt)  +  [r/T/mrm](i))J4.  We  will  assume  Tij  to  be  a  positive-definite  diagonal 
tensor,  with  Tu  <  T22.  The  isocontour  of  (r  )  is  an  ellipsoid  with  principle  axes  coincident  with 
the  axes  of  the  reference  frame,  1  /\jTaa  being  the  length  of  the  principle  axis  in  the  ea  direction. 
Thus,  Rkk(r)  is  more  correlated  in  the  x  direction  compared  to  the  y  direction  (fig  22).  Substituting 
the  correlation  into  eqn  (143)  gives 


Bij 


TijL2 


(145) 


Thus  the  eigenvalues  and  eigenvectors  of  Dij  give  the  values  and  directions  of  the  principle  axis  of 
the  isocontours  of  the  ellipsiod  characterizing  Rkk.  A  similar  conclusion  can  be  extended  for  Q*jk, 
though  we  will  not  elucidate  on  it,  since  the  nature  of  anisotropy  it  captures  is  more  complicated 
compared  to 


6.2  Formulation  of  Representation 

6.2.1  General  form  of  the  linear  representation 


In  section  6. 1  it  was  shown  that  the  structure  tensors  contain  relevant  information  about  the  anisotropy 
of  Rij.  Conversely,  we  will  assume  that  a  representation  of  the  form  FtJ(r.  B,  D,  Q*)  is  capable 

4 A  sufficient  condition  would  be  that  f(r)  is  analytical  for  all  r  €  USE,  lim^oo  f(r)  =  0.  Also,  f{r)  «  1  for 
e  <  r  <  L,  and  f{r)  <C  1  for  r  >  L  where  t<l  and  L  defi  nes  the  integral  scale  of  the  correlation. 
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D 11  <  D-22 


Figure  22:  Schematic  showing  the  connection  of  Dl3  to  the  anisotropy  of  RtJ  (r) .  An  isocontour  of 
a  hypothetical  Rkk{ r)  at  r3  =  0  is  shown  for  a  correlation  for  which  Du  <  D2 2- 


of  reproducing  the  anisotropy  of  Rij(r)  in  real  turbulent  flows.  In  order  to  ensure  that  the  model 
becomes  isotropic  when  the  structure  tensors  are  isotropic,  we  model  in  terms  of  the  normalized 
anisotropies  of  the  structure  tensors  instead.  These  anisotropic  tensors  are  given  by  b^  =  qF  — 

dij  =  'qjf-'-r  and  Q*-k  =  Thus,  we  rewrite  the  representation  as  R,j  (r)  =  q2  Frj  (r,  b,  d,  Q*), 
which  we  construct  based  on  the  following  assumptions: 


•  We  are  not  representing  the  correlation  for  the  very  small  dissipative  scales,  i.e.  we  do  not 
attempt  to  model  the  viscous  healing.  Therefore,  we  do  not  expect  FVJ  to  be  exactly  equal  to 

at  r  =  0.  Instead,  we  let  the  zero  separation  correlation  be  some  value  Brj  +  q2ABij 

•  Fij  is  assumed  to  be  linear  in  it’s  dependence  on  bij,dij  and  Q*jk.  This  is  consistent  with  the 
fact  that  the  exact  relationship  between  the  structure  tensors  and  Rij  is  linear  (eqns  (141)  - 
(143)). 

•  The  tensor  function  FtJ(r,  b,  d,  Q*)  is  invariant  to  proper  rotation  of  the  reference  frame  as 
well  as  to  changes  in  the  handedness  of  the  axes  of  the  reference  frame. 

•  The  representation  satisfies  the  symmetry  R,l3{v)  =  Rji{— r),  which  is  exactly  true  for  homo¬ 
geneous  turbulence. 


We  construct  the  most  general  linear  representation  (r,  b,  d,  Q*)  using  the  invariant  theory  of 
tensors32,  and  obtain  the  following  form  for  RtJ : 


where, 


—  Bij  +  q2 

A  Bij  +  R'j(r)  +  A(r,  b)  +  R%(  r,  d)  +  fig  W  Q) 

(146) 

Aw 

=  ./)(/•) A,  -  ,/2  (/•)/•//•, 

(147) 

4(r,b) 

=  h{r)bij  +  [fi{r)8ij  +  frl(r)rirLv  •  b  •  r 

(148) 

+h(r)[ri(r  ■  b)j  +  ry(r  •  b)4] 

A/Wd) 

=  h{r)dij  +  [/8(r)%  +  f9(r)nrj}r  •  d  •  r 

(149) 

+fio(r)[n(r  ■  d)j  +  rj{ r  •  d)<] 

AWQ*) 

=  fll{r)\^imkQklj  +  ejmkQ  kli]rlrm 

(150) 

+fi2(r)[r  j€ink  +  riejnk]Q*klmrirmrn 
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Here  fi(r)  —  f  12(f)  are  scalar  functions  of  r.  The  number  of  free  scalar  functions  /,(r)  can  be 
reduced  by  using  the  constraints  mentioned  in  section  6.2.2. 


6.2.2  Constraints 


The  Continuity  Constraint  The  continuity  equation  for  Rij(r)  in  the  context  of  homogeneous, 
incompressible  turbulence  is  given  by 


dRijjr) 

drj 


(151) 


We  get  the  following  constraints  by  substituting  eqn  (146)  into  eqn  (151),  and  equating  the  coeffi¬ 
cients  of  linearly  independent  vectors  involving  r  and  the  structure  tensors  to  zero: 


f[  +  r2f2  +  4  rf2 

=  0 

(152) 

/l  +  6  rfc  +  r2f'5  +  /' 

=  0 

(153) 

2  rfA  +  fs  +  5  rf6  +  r2  f'% 

=  0 

(154) 

fs  +  6  r/9  +  r2/'  +  f[0 

=  0 

(155) 

2r/8  +  f7  +  5r/i0  +  r2f[0 

=  0 

(156) 

r2fi2  +  7r/i2  +  f'u 

=  0 

(157) 

After  the  imposition  of  the  continuity  constraint,  R1  and  RQ  have  one  unconstrained  scalar  func¬ 
tion  left  while  Rb  and  Rd  each  have  two  unconstrained  scalar  functions  left. 


The  Self-Consistency  Constraint  When  we  evaluate  the  structure  tensors  by  substituting  the 
representation  given  by  eqn  (146)  into  the  right  hand  side  of  eqns  (147)-(150),  we  get  a  set  of 
equations  where  B,  D,  Q*  are  related  to  a  sum  of  integrals  of  the  parts  of  the  representation  given 
by  R(r  R\j,  Rfj  and  R®.  We  would  like  a  self-consistency  condition  where  when  a  structure 
tensor  is  evaluated  from  the  representation,  it  should  have  contributions  from  only  parts  of  the 
representation  relevant  to  that  structure  tensor.  More  precisely,  when  a  structure  tensor  is  zero, 
we  would  like  the  same  structure  tensor  evaluated  from  the  representation  to  be  zero.  Most  of  the 
equations  obtained  from  these  constraints  are  identically  true  due  to  the  inherent  symmetries  in 
the  representation.  The  only  equations  which  are  not  satisfied  identically  are  the  ones  obtained  by 
equating  to  zero  the  contribution  of  Rf-  to  BtJ  in  eqn  (142)  and  the  contribution  of  Rft  to  Btj  in 
eqn  (143): 


47T  7 1  ■  r  &rkdrk 

l_  r  13V4(r,b)ir 

4vr  Jv  r  d  ri  d  rj 


0 

0 


(158) 

(159) 


These  integral  constraints  are  not  strong  enough  for  the  small  separation  correlations,  because  the 
integral  does  not  have  to  be  zero  when  it  is  taken  over  finite  volumes  contained  in  V.  We  are 
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interested  in  using  these  correlations  for  LES  modeling,  and  are  therefore  interested  in  obtaining  a 
representation  of  the  correlation  for  only  small  values  of  r.  Thus,  we  apply  a  stronger  constraint 
where  we  require  that  the  angular  integral  of  the  integrands  in  eqns  (158)-(159)  are  equated  to  zero: 


1  Jft  d2Rfj(r,  d)  _ 

4tt  J  drkdrk 
1  fd2Rbkk( r,b) 

i ~J^krdZ-  = 0 


(160) 

(161) 


Where  dZ  =  sin  OdOdo  and  (r,  9,  (ft)  are  the  standard  spherical  coordinates.  Performing  the  above 
integrals  5,  we  obtain: 


2^-(15g  +  9rg'  +  r2  g"  +  35r2  f5  +  I3r3  f'5  +  r4  f£)  =  0 

^(30  f'/r  +  15/"  +  40r2/g  +  20r3/'  +  2r4/" 

+60/10  +  60r/(0  +  10r2/"  )  =  0 


(162) 


(163) 


where  g  =  3/4  +  2 /6. 


6.2.3  Power  law  model  for  /,(r) 


The  continuity  and  consistency  constraints  arise  separately  from  the  different  parts  of  the  cor¬ 
relation,  and  the  scalar  functions  are  accordingly  decoupled  from  each  other.  Examining  eqns 
(152)-(157)  and  eqns  (162)-(163)  we  see  that  it  is  possible  to  break  up  the  set  of  scalar  functions 
into  four  subsets  such  that  the  constraints  do  not  relate  functions  which  do  not  belong  to  the  same 
subset.  These  subsets  are  given  by  {/i,  f2},  {/3,  /4,  /s,  fa},  {h,  fs,  /g,  fw},  {/n,  /12}-  We  now 
postulate  that  these  functions  have  a  power  law  representation,  with  the  rationale  that  it  has  the 
ability  to  reduce  to  the  r2//3  power  law  for  isotropic  turbulence.  Also,  for  homogeneous  turbu¬ 
lence,  dependence  of  the  anisotropic  part  of  the  correlation  on  a  power  law  of  r  has  been  proposed 
theoretically31  and  reported  experimentally 33  in  the  past.  In  order  to  satisfy  the  constraints  in  a 
nontrivial  way,  the  scalar  functions  belonging  to  the  same  subset  need  to  have  the  same  power 
law,  corrected  by  an  appropriate  integral  power  of  r.  For  instance,  if  we  assume  f\  =  a^r11  and 
/2  =  a2r72,  then  eqn  (152)  yields: 


ai7i-r71  1  +  5a2r72+1  =  0  (164) 

If  71  —  1  7^  y2  +  1  then  07  =  a2  =  0  for  r  ^  0.  On  the  other  hand,  requiring  71  =  72  +  2  yields 
the  constraint  77x1  +  5a2  =  0.  Following  a  similar  procedure  for  all  the  constraints,  we  obtain 
the  form  ft(r)  =  a,rz^+n'n,  where  z[i)  are  integers,  and  oi\  —  0:4  are  the  independent  power-law 
indices  corresponding  to  each  group.  These  power-law  forms  of  fi(r)  are  then  substituted  into 
the  constraints,  given  by  eqns  (152)-(157)  and  eqns  (162)-(163)  to  yield  8  equations  for  12  af  s 
(6.4).  Therefore  there  can  be  only  four  independent  af  s,  which  we  choose  to  be  a  \ ,  a3,  a7  and  an. 

5 We  have  used  the  identities  /?•  =  ^  f  nrj  dZ.  =  and  lfjU  =  47  /  nr^n  dZ  =  +  $ikSji  + 

Siidjk)  to  derive  eqns  (162)-(163) 
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Given  a  set  of  indices  <34  —  04,  the  ratios  of  the  rest  of  the  coefficients  with  respect  to  these  four 
coefficients  are  fixed.  Thus  we  can  rewrite  eqns  (147)-(150)  as: 


ARtj(r) 


Rij(  r) 


B 


T 


=  A  Bn  +  arrai  \  r28ij  +  — r^ry 

_  ai 

Au  ,  r°4x  „2  ,  ^5  1  u  „  ,  a6.2r 
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1  a3  a3  a3 


+a7ra3{r4dji  +  [— <%r2  +  —  ryryjr  •  d  •  r  +  — r2[ry( r  ■  d)j  +  ry(r  ■  d),]} 
t  a7  07  a7  J 

Tour  ^jmkQ 

U 1 2  [fj’Gnfc  T  T'i£jnk\Qkl'mTl^mT' 


On 


(165) 


The  above  relation  implies  that  if  we  fix  (AB.  b.  d.  Q * }  then  there  are  8  degrees  of  freedom  (dof) 
in  the  representation,  given  by  ay  —  a4  and  a  1,  a3,  a7,  an. 


6.3  Fitting  the  representation  to  DNS  correlations 

6.3.1  Optimization  algorithm 


In  order  to  check  the  usefulness  of  the  representation,  we  construct  a  model  f?“odel  by  evaluating 
the  best  fit  of  the  representation  with  respect  to  the  correlation  tensor  calculated  from  DNS  of 
anisotropic  turbulent  flows,  /i'|4NS .  For  calculating  the  best  fit  we  fix  q2  and  b7l  by  using  the  value 
directly  from  DNS,  since  they  are  easily  obtained  from  the  data.  We  do  not  explicitly  vary  a7 
and  an,  and  instead  absorb  them  into  d,,:)  and  Q*jk.  Thus  the  optimal  values  of  dZJ  and  Q*jk  will 
be  correct  to  some  constant  factor.  The  search  space  for  obtaining  the  optimal  fit  is  given  by  the 
set  of  variables  (AB,  a4,  a3.  d,  Q*,p},  where  p  =  (04  +  2,  a2  +  4,  a3  +  4,  a4  +  4)  is  the  four 

dimensional  vector  of  power  law  indices.  The  best  fit  is  obtained  by  minimizing  the  L2  norm  of 

the  error  between  the  DNS  correlation  and  the  model  correlation.  This  error  is  given  by 

E{ AB,  au  a3,  d,  Q*,  a)  =  ||  ARDNS  -  ARModel||2  (166) 

where 

||T(r)||2  =  (T|T)  (167) 

and  (S|T)  =  [  .S/;( r)7';,(r)r/r  (168) 

Jv 

Sij{ r),  Ty( r)  being  tensor  fields  and  V  C  R3  having  finite  volume.  V'  needs  to  have  finite 
volume  since  our  model  for  Rij(r)  does  not  decay  to  zero  for  large  r  and  therefore  the  error 
norm  E  depends  implicitly  on  our  choice  of  the  volume  V'.  We  will  however  perform  our  error 
minimization  over  a  fixed  V ,  which  will  be  chosen  over  an  appropriate  inertial  range.  For  per¬ 
forming  the  minimization,  we  take  advantage  of  the  fact  that  the  dependence  of  ARModel  on 
(AB,  ai,  a3,  d,  Q*}  is  linear,  which  allows  to  use  a  two  step  approach.  In  the  first  step,  we  use 
least  squares  minimization  to  obtain  the  minimum  E  for  a  given  p,  i.e.  we  obtain  E'{ p7)  = 
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E( AB,  a?,  ag,  d°,  Q*°,  p')  such  that  for  a  fixed  p',  the  values  AB°,  a^a^d0  and  Q*°  are  the 
solutions  to  dE/dABij  =  0,  dE/dai  =  0,  d E/da3  =  0,  dE/ddij  =  0  and  dE / dQ*jk  =  0 
respectively.  If  the  flow  does  not  contain  any  symmetries  (other  than  homogeniety),  then  this 
step  involves  optimization  over  20  variables  6.  In  the  second  step,  we  minimize  E'( p)  over 
(p;  Pi  £  [— 1  +  c0,  oo],  i  =  1,  2,  3, 4}7,  where  c0  is  a  small  positive  constant,  chosen  to  be  0.03 
in  our  case.  The  second  step  is  carried  out  using  the  PRAXIS  optimization  software  package, 
which  is  based  on  ’Brent’s  method’ 34 .  If  p  =  pmin  is  the  exact  minima  for  ^'(p),  then  the  settings 
in  the  PRAXIS  software  ensure  that  the  optimal  p  is  calculated  within  an  error  e  =  0.01  of  pmin. 


6.3.2  Results  for  homogeneous  turbulence 


The  first  fit  of  our  representation  is  done  with  respect  to  correlations  constructed  from  the  DNS 
velocity  field  of  an  initially  isotropic  and  homogeneous  turbulent  flow  subjected  to  mean  shear 
(uij)  =  8ii8j2S.  The  details  of  the  DNS  can  be  found  in35.  Some  of  the  relevant  details  of  the  field 
used  are  given  in  table  5.  The  volume  over  which  the  fit  is  performed  is  given  by  the  spherically 
symmetric  annulus  V'  =  (r;  0.26  <  r  <  0.6},  where  the  lower  and  upper  limit  of  r  correspond  to 
the  Taylor  micro-scale  and  the  integral  length  scale  respectively.  Due  to  the  symmetries  in  the  flow, 
several  components  of  (AB,d,  Q*}  are  identically  zero,  resulting,  for  a  fixed  p,  in  12  degrees 
of  freedom  in  the  representation.  The  relative  error  obtained  between  the  best  fit  and  the  DNS 
correlation,  given  by  ||ARModel  —  ARDNS||/||  ARDNS||  is  10%.  A  more  qualitative  comparison  is 
done  in  figure  23,  where  we  have  chosen  to  show  the  isocontours  of  r3  =  0  because  the  correlations 
have  the  least  symmetry  in  this  plane. 

We  take  note  of  two  major  differences  between  the  model  and  the  DNS.  Firstly,  the  isocontours 
of  AR22  in  the  model  have  a  different  shape  when  compared  with  that  of  the  DNS.  This  occurs 
because  the  Reynolds  stress  i?2 2(0)  is  much  smaller  than  the  other  normal  Reynolds  stress  compo¬ 
nents  (table  5),  and  therefore  this  component  does  not  get  as  good  a  fit  as  the  other  components. 
Secondly,  we  notice  that  the  isocontours  of  /?[}Tlodel  are  not  as  elongated  as  the  isocontours  of  f?°NS 
in  the  streamwise  direction.  The  elongation  in  /?II)1NS  happens  due  to  the  fact  that  the  high  mean 
shear  picks  out  the  streamwise  axis  of  the  isocontours  of  the  Rn  component  to  grow  at  a  faster  rate 
compared  to  other  components36.  However,  the  volume  V'  over  which  the  best  fit  is  conducted  is 
spherically  symmetric,  and  therefore  tries  to  fit  the  same  power  law  functions  upto  an  equal  extent 
in  all  directions.  This  causes  the  fall-off  in  the  correlation  to  get  averaged  over  all  directions  and 
components,  and  thus  the  excessive  elongation  in  the  streamwise  direction  for  Rn  does  not  get 
properly  captured. 

On  the  other  hand,  we  also  observe  that  for  r  e  V’,  the  isocontours  of  the  model  and  DNS  cor¬ 
relations  seem  to  have  their  principle  axes  inclined  by  the  same  amounts  for  all  the  components 
(i.e.  except  R2 2).  Also,  specific  details  of  the  isocontours  of  each  component  (e.g.  the  elliptic 

6The  20  variables  arise  after  taking  into  account  the  symmetries  of  the  tensors  in  the  search  space.  There  are  6 
independent  components  in  which  has  the  symmetry  A Bij  =  A Bji,  1  component  each  from  a\  and  013,  5 

independent  components  from  dtl  which  has  the  symmetries  j  =  dji  and  dkk  =  0,  7  independent  components  from 
Q*jk,  which  is  invariant  to  any  permutations  in  it’s  indices  and  has  the  symmetry  Q*kk  =  0 

7The  lower  limit  of  p  is  needed  to  ensure  that  while  performing  least  square  error  minimization,  inner 
products  amongst  tensor  components  of  the  form  given  in  eqn  (168)  do  not  diverge.  That  is,  Vr  £  V 
r2  f  Sij(r)Sij(r)dZ  =  0(1)  and  r2  Ty  (r)Tij(r)dZ  =  0(1)  =>  r2  §  Sij(r)Tij(r)dZ  =  0(1)  (from  the  fact 
that  (S,  T)  =  Sij (r)Tij (r)dZ  is  also  an  inner  product  and  the  Cauchy-Schwartz  inequality) 
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Table  5:  Specifics  of  the  homogeneous  turbulent  field  (MKS  units) 


Resolution 

2563 

Box  Size  (lx  x  ly  x  lz) 

2tt  x  2tt  x  2n 

Integral  length  scale  L  (3n  x  f  k^1E(k)dk/  f  E(k)dk) 

0.6 

Taylor  length  scale  A 

0.26 

rms  velocity  (q) 

(W«i) ,  («i« i) ,  W2U2)  >  W3U3)) 

6.37 

(26.12,-5.69,5.55,8.92) 

Re  a  at  St  —  0  (i.e.  of  initially  isotropic  field) 

45.35 

Sq2/e  at  St  =  0 

35.85 

isocontours  for  the  normal  components  and  the  dual-lobed  shape  of  the  isocontours  of  f?i2(r))  are 
captured  well  by  the  model. 


Ideally,  we  would  also  like  to  compare  the  values  of  the  structure  functions  calculated  from  the 
DNS  and  Model  correlations  using  eqns  (141)-(143).  However,  the  model  assumes  a  power  law 
dependence  for  r,  and  is  therefore  not  valid  for  all  r  e  M3  (specifically  for  r  much  smaller  than  the 
Taylor  micro-scale  and  much  larger  than  the  integral  scales).  Instead,  we  define  a  scale-dependent 
dimensionality  tensor  D'VJ(r\h  r )  given  by: 


D'ij(ro,r) 


2_  f  1  d2Rkk(iJ)djJ 
4vr  Jv"(r0,r)  r'  dr'fir'j 


(169) 


Where  V"(ro,  r)  =  (r';  r0  <  r'  <  r}.  From  figure  24  we  see  that  there  is  excellent  agreement  for 
D\ A,  r)  for  A  <  r  <  L.  Thus,  our  model,  supplemented  with  a  good  model  for  R,3  at  very  large 
and  very  small  separations  will  ensure  that  we  can  obtain  correct  values  of  DV}  from  the  correlation. 


6.3.3  Results  for  turbulent  channel  flow 


The  second  set  of  correlations  used  to  fit  our  representation  were  calculated  from  DNS  fields  of 
fully  developed  turbulent  channel  flow  at  Rer  =  940.  Unlike  the  previous  case  this  is  a  highly 
inhomogenous  flow  and  we  try  to  examine  the  extent  to  which  our  homogenous  representation  will 
be  valid.  The  main  reason  to  carry  out  this  exercise  is  to  find  out  whether  the  tensors  forms  in  eqns 
(147)-(150)  need  to  depend  on  any  additional  tensors  due  to  the  presence  of  inhomogeniety. 

The  only  direction  of  inhomogeniety  is  the  one  perpendicular  to  the  channel  wall.  The  x  axis  is 
aligned  in  the  streamwise  direction  and  the  y  axis  is  aligned  in  the  direction  of  inhomogeniety. 
Therefore  the  correlation  will  be  redefined  as  Rij(y,  r)  =  (u'(x)u' (x  +  r))  | X2=y.  Our  model  does 
not  satisfy  the  continuity  equation  8  or  the  symmetry  of  the  correlation9  exactly,  but  for  small  r/y 
values  we  expect  the  effects  of  inhomogeniety  to  be  small  and  we  will  assume  that  Rij(y,  r)  is 
locally  homogenous  at  y.  This  of  course  implies  that  all  our  fitting  parameters  are  functions  of  y. 

After  performing  the  best  fit  optimization,  we  found  the  relative  error  between  our  model  and  the 
fit  to  be  6%.  Figure  25  shows  that  the  model  fits  the  data  quite  well  for  the  reference  y+  =  220 

8Given  by  (^Si2  -  ^r)  Rij(y,  r)  =  £-Ri3(y,r)  =  0 

9RtJ(y,r)  =  Rji(y  +  r2,-r) 
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Figure  23:  Isocontours  of  Ai2y(r)  =  —  Rrj{()),  plotted  for  r2  =  0.  The  correlations  in 

the  left  column  were  calculated  from  DNS  of  homogeneous  shear  flow,  and  the  correlation  in  the 
right  column  from  the  best  fit  of  the  model  with  DNS.  The  contour  levels  for  given  a  correlation 
component  (i.e.  for  the  same  row)  have  the  same  range. 
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Figure  24:  Comparison  of  D'?(r0,  r)  —  ^  f  where  V"(r0,  r )  =  {r';  r0  <  r'  < 

V"(ro,r)  *  7 

r}.  Here  r0  =  A  =  0.26,  and  has  calculated  from  both  model  (dashed  lines)  and  DNS  (bold 
lines). 


location  in  the  log  layer.  The  inclination  of  principle  axes  are  well  captured  and  the  magnitude  of 
the  contour  levels  match.  The  main  shortcoming  of  our  model  lies  in  it’s  inability  to  capture  the 
effects  of  inhomogeniety,  which  is  best  seen  in  the  R12  component,  where  the  isocontours  of  the 
DNS  extend  out  further  in  the  positive  r2  half  of  the  plane  as  compared  to  the  negative  r2  half  due  to 
the  presence  of  the  wall.  Our  model  on  the  other  hand  is  symmetric  in  r,  and  is  not  able  to  capture 
this  aspect  of  the  data,  which  points  to  the  need  for  further  improvements  in  the  representation. 


6.4  Conclusions 

A  finite  dimensional  representation  of  the  two  point  correlation  tensor  was  constructed  using  the 
invariant  theory  of  tensors.  A  novel  approach  was  taken  in  this  work,  where  the  representation 
depended  upon  structure  tensors.  Despite  being  single-point  moments,  these  tensors  capture  the 
extent  and  orientation  of  the  elongation  of  Rkk  for  values  of  r  of  the  order  of  the  integral  scale  of 
the  correlation.  Our  representation  depended  linearly  on  the  structure  tensors,  and  we  were  able  to 
significantly  reduce  the  number  of  degrees  of  freedom  in  our  model  using  the  continuity  and  self 
consistency  constraints.  A  power-law  dependence  on  r  was  assumed  in  the  model,  which  prevented 
us  from  modeling  the  dissipative  range  and  the  large-separation  range.  In  order  to  fit  the  correlation 
to  the  DNS  data,  we  had  to  allow  the  Reynolds  stress  of  the  model  to  be  a  degree  of  freedom 
for  our  fit.  The  results  show  that  our  model  fits  the  correlations  calculated  from  DNS  data  well 
both  quantitatively  and  qualitatively,  and  that  complemented  with  a  model  for  the  large-separation 
and  dissipative  range  separation,  one  could  extract  the  structure  tensors  from  experimental  and 
numerical  measurements  of  Rij  by  fitting  it  to  our  model.  In  fact  this  model  enables  us  to  construct 
good  reconstruction  of  R^  using  a  finite  number  of  data  points,  which  can  be  potentially  useful  for 
OLES  formulations,  where  information  about  the  correlation  is  required  from  a  finite  dimensional 
field.  However,  the  model  does  require  a  variety  of  improvements,  most  urgent  being  the  need  to 
model  dissipative  range  and  large-separation  range  of  R^,  and  constructing  a  representation  for 
the  anisotropy  which  is  valid  for  inhomogenous  flows. 
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Figure  25:  Isocontours  of  A Rij(y,  r)  =  Rij(y,  r)  —  Rij(y,  0)  at  y+  =  ReTy/h  =  220,  plotted  for 
rz  =  0.  The  correlations  in  the  left  column  were  calculated  from  DNS  of  channel  flow,  and  the 
correlation  in  the  right  column  from  the  best  fit  of  the  model  with  DNS.  The  contour  levels  for 
given  a  correlation  component  (i.e.  for  the  same  row)  have  the  same  range. 
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Constraints  arising  from  /.t(r)  «  rz\A+an 

In  section  6.2.3  it  was  observed  that  the  set  of  functions  {fi(r)}  can  be  divided  into  four  disjoint 
subsets,  according  the  the  set  of  constraints  satisfied  by  them.  These  subsets  are  given  by  {/ 1,  /2}, 
{hi  h,  h,  fe},  {/7,  hi  hi  fio},  {/11,  /12}.  Assuming  a  power  law  form  /;(r)  =  a.;r7i,  it  is  then 
necessary  for  7*’s  to  have  a  form  so  that  every  term  in  a  given  constraint  equation  has  the  same 
power  of  r  (otherwise  we  get  the  solution  a*  =  0  for  all  V s).  We  choose  /,’s  to  have  the  form: 


h 

=  air2+ai 

f2  =  a2rai 

(170) 

h 

=  a3r4+“2 

U  =  «4r2+a2 

=  abr°2 

h- 

=  a6r2+a2 

(171) 

h 

=  a7r4+"3 

h  =  a8r2+Q3 

h 

=  a9r“3 

fio 

=  ai0r2+Q3 

(172) 

fu 

=  an-r2+“4 

fi2  =  ai2r°4 

(173) 

Substituting  the  forms  of  f\  —  /12  in  eqns  ( 170)-(  173)  into  eqns  (152)-(157)  and  eqns  (162)-(163) 
we  get  the  following  set  of  linear  equations  for  a\  —  a[2: 

ai  (2  +  01)  +  a2(4  +  01)  =  0  (174) 

04(2  +  o2)  +  a6(2  +  o2)  +  05(6  +  a2)  =  0 

2a4  +  03(4  +  o2)  +  ag(7  +  o2)  =  0 

3a4  T  a 5  T  2cig  =  0 

a8(2  +  0:3)  +  ai0(2  +  0:3)  +  ag(6  +  0:3)  =  0 

2a8  +  07(4  +  0:3)  +  040(7  +  0:3)  =  0 

(15a7  +  2(ag  +  5aio))(4  +  0:3)  (5  +  0:3)  =  0 

On(2  +  0:4)  +  Oi2(7  +  0:4)  =  0 

Choosing  04,  a3,  07,  an  to  be  the  independent  coefficients,  the  ratios  of  the  remaining  coefficients 
with  respect  to  the  four  independent  coefficients  become  fixed,  and  are  given  as  follows: 

a2 

Ol 

a  5 
03 
ci  8 

CLy 

a  10 

Cij 

As  long  as  the  denominators  of  eqns  (175)-(178)  do  not  vanish10,  the  ratios  of  coefficients  are 
uniquely  determined.  From  section  6.3.1  we  know  that  the  search  space  for  p  =  (07  +  2,  a2  + 
4,  o3  +  4,  o4  +  4)  is  given  by  (p;  pi  e  [—1  +  c0,  00],  i  =  1,  2,  3, 4},  where  c0  =  0.03.  Therefore 
none  of  the  roots  of  the  denominators  of  eqns  (175)-(178)  are  within  the  range  of  the  search  space. 

1()This  is  true  if  04  ^  —4,  a2  ^  —7  ±  \/3,  <43  ^  —10,  a 4  ^  —7 


(2  +  Q4)  04  _  1  / 

(4  +  cti)  (is  2  \ 

(2  +  ct2)(4  +  a2) 

2(46  +  ct2(14  +  a2))  ’ 

(l _ *_) 

(10  +  a3) ) 

_  (14  +  CC3)  a12 
(10  +  0:3)’  an 


V  (46  +  o2(14  +  o2))y 
Qjq  (4  +  a2)(8  +  o2) 

a3  (46  +  a2(14  +  o2)) 

a 9  _  —5(2  +  03) 

cij  2(10  +  03) 

(2  +  04) 

(7  +  04) 


(175) 


(176) 


(177) 


(178) 
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7  Testing  and  Evaluation  of  the  Filtered  Wall  Formulation  for 
the  Wall-Layer 


One  of  the  pacing  problems  in  the  development  of  reliable  large  eddy  simulation  (LES)  models 
for  use  in  turbulent  flows  of  technological  interest  is  the  so-called  LES  wall-modeling  problem  37. 
It  arises  because  the  length-scale  associated  with  the  wall  layer  of  a  turbulent  wall-bounded  shear 
flow  (wall  units)  gets  smaller  relative  to  the  shear  layer  thickness  approximately  like  the  inverse 
Reynolds  number  (like  Rer  7^8  in  the  channel  flow).  The  “large-scale”  turbulence  in  this  thin 
layer  also  scales  in  wall  units.  If  the  cost  of  an  LES  of  wall-bounded  flows  is  to  remain  finite 
in  the  limit  of  infinite  Reynolds  number,  then  this  wall  layer  and  the  large-scale  turbulence  it 
supports  cannot  be  represented  directly,  and  so  must  be  modeled.  However,  current  LES  models 
are  generally  not  valid  for  this  near- wall  layer  because  underlying  assumptions  such  as  small-scale 
homogeneity  and  isotropy  are  not  valid.  The  alternative  is  to  resolve  the  near- wall  turbulence.  The 
most  successful  LES  of  wall-bounded  shear  flows  employ  this  technique,  though  this  is  clearly  not 
viable  for  arbitrarily  large  Reynolds  number.  In  this  paper,  we  propose  a  possible  solution  to  this 
wall-modeling  problem. 

Our  approach  is  motivated  by  the  observation  that  in  an  LES,  locating  anything,  including  the 
wall,  to  more  precision  than  the  filter  width  is  inconsistent  with  the  representation.  This  leads 
us  to  a  formulation  in  which  the  wall  is  filtered  as  well  as  the  turbulence.  This  allows  the  use  of 
homogeneous  or  nearly  homogeneous  filters  normal  to  the  wall,  avoiding  commutation  error.  Also, 
the  the  dynamics  of  the  sublayer  can  be  effectively  filtered  out,  which  we  will  see  is  advantageous. 
Filtering  the  wall,  however,  also  introduces  an  extra  term  in  the  filtered  equations  which  must  be 
modeled,  along  with  the  usual  subgrid  stress  term.  A  modeling  approach  for  this  extra  term  is 
proposed  here. 

To  allow  the  filtered  boundary  approach  to  be  evaluated  with  minimum  uncertainty  arising  from  the 
subgrid  stress  model,  optimal  LES  models  will  be  used1,8’26-28.  In  optimal  LES,  the  subgrid  force 
term  (or  the  subgrid  stress)  is  approximated  using  stochastic  estimation.  Optimal  LES  is  a  formal 
approximation  to  what  we  have  called  the  ideal  LES  evolution  \  which  can  be  shown  to  produce 
one-time  statistics  that  are  exact,  and  minimum  mean-square  variation  in  the  instantaneous  large- 
scale  evolution.  The  optimal  LES  formalism  has  the  advantage  in  this  context  of  being  valid  even 
in  the  absence  of  small-scale  isotropy  or  homogeneity;  that  is,  it  is  valid  for  near-wall  turbulence. 
As  input,  optimal  LES  requires  detailed  two-point  correlation  data.  For  the  purposes  of  testing 
the  viability  of  the  proposed  wall-modeling  approach,  this  data  has  been  obtained  from  the  direct 
numerical  simulation  data  of  Moser  et  al 38 . 

In  the  remainder  of  this  paper,  the  filtered  boundary  formulation  is  introduced  and  a  test  of  it’s 
capabilities  is  presented.  The  optimal  LES  models  used  here  are  briefly  described.  The  results  of 
filtered  boundary  LES  of  the  turbulent  channel  at  ReT  =  590  are  then  presented  followed  by  a 
brief  discussion  of  the  implications  of  these  results,  particularly  the  properties  of  the  subgrid  stress 
model  needed  in  this  formulation. 
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7.1  Filtered  Wall  Formulation 


In  the  filtered  boundary  LES  formulation,  the  wall-bounded  domain  is  embedded  in  a  larger  do¬ 
main,  with  the  Navier-Stokes  equations  applied  to  the  interior,  and  u  =  0  applied  to  the  exterior 
domains.  However,  since  u  =  0  is  a  solution  to  the  Navier-Stokes  equations,  we  can  also  extend 
Navier-Stokes  into  the  exterior  domain,  with  the  wall  boundary  supporting  a  stress  discontinuity 
just  sufficient  to  maintain  u  =  0  in  the  exterior.  A  filter  is  then  applied  to  this  larger  domain.  The 
resulting  filtered  equations  are  then: 


duj 
dt  + 


dujUi 

dxj 


dp  1  dui 

dxi  Re  dxjdxj 


+  b{  +  Mi 


(179) 


where  Mi  =  —drij/dxi  is  the  usual  LES  model  term  and  b,  is  the  boundary  term.  The  subgrid 
stress  is  given  by  tV)  =  iRuj  —  UiUj  and  the  boundary  term  (b,)  can  be  written 


bi(x)  —  /  <Jij(xr)rijG(x  —  x')  dx 

JdR 

where  a  is  the  stress  at  the  boundary,  including  pressure  and  viscous  stress,  dR  is  the  boundary  of 
the  fluid  region  R  and  rij  is  the  unit  normal  to  the  surface.  If  a  sufficiently  fine  filter  width  is  used 
then  Mi  is  negligible  and  the  only  effect  is  the  filtering  of  the  boundary  (i.e.  a  filtered  boundary 
“DNS”).  Such  a  “DNS”  was  used  as  a  test  case  (see  below).  In  many  LES  of  wall  bounded 
flows,  approximate  boundary  conditions  are  used  to  model  the  effect  of  the  wall  layer39.  The 
approximate  boundary  conditions  are  prescribed  in  terms  of  the  wall  shear  stress,  so  wall  stresses 
must  be  determined  in  terms  of  the  resolved  velocities.  In  the  present  formulation,  the  unfiltered 
wall  stresses  are  also  required,  and  for  analogues  reason. 


In  the  current  description,  in  which  the  unfiltered  velocity  is  zero  in  the  buffer  domain  external 
to  the  walls,  the  wall  stress  is  the  surface  forcing  required  to  ensure  that  momentum  and  energy 
are  not  transfered  to  the  buffer  domain.  That  is,  that  the  velocity  remains  zero.  This  suggests  a 
technique  for  determining  the  wall  stress.  It  would  be  inappropriate  to  define  a  force  to  make  the 
velocity  zero  at  the  boundary  as  in  embedded  boundary  numerical  methods4041,  since  the  filtered 
velocity  is  not  zero  at  the  wall.  Instead,  we  choose  awau  to  minimize  the  transport  of  momentum 
to  the  exterior  domain.  To  this  end,  the  wall  stresses  at  each  time  step  are  defined  by  minimizing 

E  —  /  |u|J  +  a 

Jb 


Ou 

~dt 


dx 


(180) 


where  the  integral  is  over  the  buffer  domain.  The  |u|2  term  forces  the  energy  in  the  buffer  domain 


to  be  small,  and  the  a 
taneously  correct  smal 


Ou 

W 


term  ensures  that  large  errors  are  not  incurred  in  the  stress  to  instan- 
errors  in  the  buffer  domain.  The  constant  a  controls  the  balance  between 


these  two  competing  requirements  and  is  set  to  a  value  of  order  At2. 


In  this  paper,  we  consider  turbulent  flow  in  a  channel  and  apply  a  Fourier  cut-off  filter  and  Fourier 
spectral  method  in  all  three  spatial  directions.  In  this  case,  the  minimization  of  E  required  to 
determine  the  wall  stress  a  is  straight  forward  since  it  can  be  done  independently  for  each  (kx,  kz) 
wavenumber,  resulting  in  a  6-parameter  optimization  in  ( crxy ,  ayy,  crZy)42- 
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Figure  26:  Effect  of  the  boundary  terms  in  the  evolution  of  small  disturbances  in  a  channel  flow,  (a) 
filtered  u  velocity,  (b)  filtered  v  velocity,  (c)  filtered  pressure,  (d)  pressure  gradient,  (e)  boundary 
term  for  v  equation,  (f)  pressure  gradient  +  boundary  term.  —  real  part, - imaginary  part. 


To  evaluate  this  approach,  we  consider  two  test  cases:  propagation  of  an  Orr-Sommerfeld  wave 
and  low  Reynolds  number  turbulence  in  a  channel.  In  both  cases,  the  Fourier  cut-off  filter  is  fine 
enough  to  make  the  model  term  Mi  negligible. 

In  the  Orr-Sommerfeld  case,  the  simulated  growth  rate  was  within  0.25%  of  the  exact  value  for 
the  case  considered.  More  interesting  is  the  role  the  boundary  term  plays.  Consider  the  exact 
unfiltered  pressure  fluctuations.  They  are  formally  zero  in  the  exterior,  resulting  in  a  discontinuity 
in  pressure,  and  the  resulting  Gibbs  phenomenon  in  the  filtered  pressure  is  shown  in  figure  26c.  The 
wall  normal  pressure  gradient  appears  in  the  (-momentum  equation,  and  this  quantity  is  dominated 
by  the  filtered  delta  function  at  the  boundary  and  the  resulting  Gibbs  phenomenon  (figure  26d).  Yet 
the  Gibbs  phenomenon  in  velocity  perturbations  in  figure  26a  and  b  is  imperceptible.  The  reason 
is  that  the  term  bv  (figure  26e)  has  exactly  the  same  structure  as  the  pressure  gradient  and  cancels 
the  Gibbs  phenomenon  (figure  26f).  The  role  of  the  boundary  terms  in  the  momentum  equation  is 
thus  to  regularize  the  stress  discontinuities  at  the  wall  (both  pressure  and  viscous  stresses). 

To  assess  the  applicability  of  this  technique  in  simulating  turbulent  flow,  a  fully  developed  channel 
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20 , 


Figure  27:  Mean  (left)  and  rms  (right)  velocities  in  turbulent  channel  flow  at  Rer  =  180. 
- present  using  filtered  boundary  formulation, - Moser  et.  al.  (1999) 


flow  is  computed  on  a  128  x  256  x  128  grid  in  x,  y  and  z  respectively,  with  20  y-points  in  the  buffer 
region.  The  friction  Reynolds  number  is  ReT  =  180  and  the  domain  size  is  the  same  as  in  the  DNS 
of  Moser  et  at3S.  The  mean  and  rms  velocities  from  this  simulation  are  in  excellent  agreement  with 
those  from  the  DNS  (see  figure  27),  and  the  near  wall  turbulence  exhibits  the  familiar  structures, 
such  as  streaks  and  inclined  shear  layers. 


7.2  Optimal  Large  Eddy  Simulation 

Optimal  LES  is  based  on  the  observation  that  there  is  an  ideal  LES  model,  which  guarantees  correct 
single  time  statistics  and  minimum  error  is  short-time  dynamics 1,23 .  This  ideal  model  given  by 

mi={Mi\u  =  w),  (181) 

where  rnt  is  the  model  for  the  term  M%  in  (179),  w  is  the  LES  field,  and  u  is  the  filtered  real 
turbulence.  In  essence,  this  is  the  average  of  Mi  over  all  turbulence  fields  that  map  to  the  LES 
field  through  the  filter.  Unfortunately  this  model  is  intractable,  so  in  optimal  LES  we  approximate 
this  model  using  stochastic  estimation 3,443 .  In  the  LES  performed  here,  the  stochastic  estimation 
formulation  is  simplified  by  the  homogeneity  of  the  channel  flow  in  directions  parallel  to  the  wall, 
and  the  formulation  must  be  further  simplified  to  avoid  problems  of  over  generalization28.  The 
linear  stochastic  estimate  used  here  can  thus  be  written: 

A  i(y)  =  (Mi)  +  kij(y)Ej(y)  (182) 

(Ml(y)Et(y))  =  kij(y){Ej(y)El(y))  (183) 

where  :  indicates  the  Fourier  transform,  and  the  event  vector  E:I  is  a  vector  consisting  of  the  fluc¬ 
tuating  LES  velocities  w'  and  their  y  derivatives.  The  correlations  appearing  in  (183)  must  be 
determined  to  complete  the  model.  For  the  purposes  of  the  test  described  below,  the  correlations 
were  evaluated  using  the  DNS  data  from  Moser  et  al3S  at  ReT  =  590.  Using  DNS  data  allows  the 
optimal  LES  formulation  to  be  evaluated  without  uncertainties  introduced  by  further  modeling  of 
the  correlations,  or  other  modeling  approximations. 
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Figure  28:  Mean  velocity  (left)  and  rms  streamwise  velocity  fluctuations  (right)  in  a  turbulent 
channel  at  ReT  =  590.  The  LES  was  performed  using  the  filtered  boundary  optimal  LES  formula¬ 
tion. 

7.3  Filtered- Wall  LES  Results 

The  filtered  boundary  formulation  and  the  optimal  LES  model  were  used  to  perform  an  LES  of 
turbulent  channel  flow  with  bulk  Reynolds  number  Reb  =  10,950  corresponding  to  a  channel 
with  ReT  =  590.  Periodic  boundary  conditions  were  used  in  streamwise  (x)  and  spanwise  (z) 
directions,  with  domain  sizes  Lx  =  2nh  and  Lz  =  irh  (h  is  the  channel  half  width).  DNS  of  this 
case  was  performed  by  Moser  et  al3S,  and  optimal  LES  were  performed  by  Volker  et  al 28. 

To  accommodate  the  filtered  boundary  formulation,  a  buffer  region  is  added  outside  the  channel 
and  periodic  boundary  conditions  are  used  in  the  extended  wall-normal  ( y )  domain.  Fourier  cut-off 
filters  with  effective  filter  widths  of  Ax+  =  116,  A y+  =  37  and  A z+  =  58  are  used  in  the  three 
spatial  directions.  In  x  and  z,  these  are  the  same  filters  used  in  Volker  et  al 28 .  Note  that  these  filter 
widths  are  sufficiently  large  to  eliminate  the  structure  of  the  near-wall  viscous  and  buffer  layers. 

The  filtered  boundary  model  and  the  optimal  LES  model  were  used  to  perform  an  LES  of  the 
channel  flow.  The  statistical  correlations  required  as  input  to  the  optimal  LES  formulation  were 
determined  from  the  DNS  of  Moser  et  al 38.  Sample  results  from  this  simulation  are  shown  in 
figure  28.  Note  that  despite  the  fact  that  the  wall  layer  was  not  resolved,  both  the  mean  velocity 
and  the  rms  velocities  are  in  remarkably  good  agreement  with  the  filtered  DNS.  Furthermore,  the 
one-dimensional  spectra  in  the  x  and  z  directions  are  in  good  agreement  with  the  filtered  DNS 
data. 

It  should  be  noted  that  the  optimal  LES  form  used  here  is  the  simplest  of  those  proposed  by  Volker 
et  al 28  for  the  channel,  and  that  in  that  study  models  of  this  form  performed  poorly.  The  reason  was 
that  this  form  does  not  properly  represent  the  wall-normal  transport  of  energy  and  Reynolds  stress. 
As  pointed  out  by  Haertel  et  al 44,  in  the  absence  of  wall-normal  filtering,  the  contribution  of  the 
subgrid  term  to  the  resolved-scale  energy  equation  is  positive  near  the  wall  (see  figure  29),  which 
is  due  to  the  subgrid  contribution  to  the  transport  of  energy  from  the  production  peak  to  toward 
the  wall.  Volker  et  al 28  found  that  a  more  complicated  form  that  did  represent  the  wall  normal 
transport  produced  a  model  that  performed  very  well. 

However,  with  coarse  wall-normal  filtering,  the  subgrid  term  is  everywhere  dissipative,  as  is  also 
shown  in  figure  29.  There  is  therefore  no  instability  introduced  by  using  the  simple  model  form 
used  here,  with  the  result  that  the  model  produces  high  quality  simulations,  as  shown  above. 
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Figure  29:  Subgrid  energy  transfer  (■ )  in  turbulent  channel  flow,  as  measured  from  DNS,38 
with  and  without  wall-normal  filtering,  and  a  posteriori  from  the  filtered- wall  LES. 

Also  notice  the  large  fluctuations  (wiggles)  in  the  profile  of  the  subgrid  dissipation  shown  in  fig¬ 
ure  29.  The  a  posteriori  dissipation  from  the  LES  reproduces  these  features,  and  is  otherwise  in 
general  agreement  with  the  a  priori  values.  The  reason  for  these  oscillations  is  the  Gibbs  phe¬ 
nomenon  introduced  by  filtering  through  the  wall,  and  sublayer.  These  fluctuations  are  apparently 
important  in  the  simulation,  because  they  cancel  similar  Gibbs  phenomena  appearing  in  other 
terms,  as  with  the  Orr-Sommerfeld  test  case  described  in  section  2. 


7.4  Discussion 

The  results  described  above  are  intriguing  because  they  suggest  that  it  is  not  necessary  to  resolve 
the  near- wall  layer  in  an  LES  to  obtain  an  accurate  simulation  of  a  wall-bounded  flow.  However, 
because  the  simulations  reported  here  were  based  on  knowledge  of  statistical  correlations  obtained 
from  DNS,  the  work  presented  here  does  not  constitute  a  practical  broadly-applicable  LES  model. 
For  this,  the  need  for  DNS  statistical  data  must  be  overcome.  None-the-less,  the  current  results  do 
demonstrate  the  value  of  the  wall-filtering  approach,  and  the  “no-leakage”  optimization  model  for 
the  wall  stresses.  It  would  appear  that  this  approach  may  form  the  foundation  of  a  solution  to  the 
well-known  LES  wall  modeling  problem. 

The  results  of  this  study  are  also  curious,  because  the  “no-leakage”  model  does  not  incorporate  any 
information  regarding  the  properties  of  the  near-wall  sublayer  that  has  been  filtered  away.  Whereas 
it  has  been  commonly  assumed  that  the  wall  boundary  treatment  must  account  for  the  presence  of 
the  sublayer.  Further,  the  “no-leakage”  model  cannot  by  its  nature  account  for  the  presence  or 
absence  of  roughness.  This  raises  the  question  of  how  information  regarding  near-wall  properties 
of  the  turbulence  (which  are  directly  affected  by  roughness)  is  incorporated  into  the  LES.  The 
answer  is  that  the  nature  of  the  near-wall  turbulence  is  reflected  in  the  subgrid  fluctuations,  which 
affect  the  LES  through  the  sub  grid  model  term  M,  =  —drij/dxj.  In  the  simulation  reported  here, 
the  optimal  LES  model  of  the  subgrid  term  is  formulated  based  on  statistical  data  obtained  from  a 
smooth- wall  DNS.  Thus,  the  simulation  produced  results  consistent  with  a  smooth- walled  channel. 
We  propose  that  an  optimal  model  formulated  for  a  particular  wall  roughness  would  also  produce 
an  LES  consistent  with  that  rough  wall. 

It  appears  then,  that  the  LES  wall  boundary  condition  problem  may  not  be  as  difficult  as  has 
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been  commonly  believed.  The  filtered  boundary,  “no-leakage”  formulation,  appears  to  provide  an 
adequate  boundary  condition.  Instead,  the  challenge  is  the  modeling  of  the  volumetric  subgrid 
stress  term.  However,  at  high  Reynolds  number,  and  with  the  sublayer  effectively  filtered  out, 
the  model  need  only  represent  the  subgrid  terms  in  the  log-layer,  where  similarity  scaling  makes 
the  modeling  easier.  As  indicated  above,  it  is  also  necessary  that  the  effect  of  roughness  on  the 
log-layer  subgrid  model  be  characterized. 

To  relieve  the  need  for  DNS  data,  which  was  used  for  subgrid  modeling  here,  a  model  for  the 
subgrid  term  with  minimal  empirical  input  and  which  retains  the  good  properties  of  the  optimal 
LES  model  is  needed.  But,  what  properties  of  the  subgrid  model  used  here  are  important  for  LES 
performance  in  wall-bounded  flows?  The  formulation  of  the  optimal  model  provides  guidance, 
because  it  was  designed  to  accurately  represent  only  the  transfer  of  energy  and  Reynolds  stress 
to  the  small  scales,  and  the  direct  contribution  of  the  subgrid  to  the  mean  Reynolds  stress.  More 
precisely,  it  provides  an  accurate  representation  of  the  dependence  of  the  subgrid  dissipation  on 
the  x  and  z  wavenumbers  and  y  location,  as  well  as  the  anisotropy  of  the  dissipation.  Whether  this 
much  detail  in  the  representation  of  the  subgrid  dissipation  is  needed,  is  not  at  this  point  clear,  and 
is  being  investigated.  Perhaps,  if  less  detail  is  sufficient,  standard  models  (e.g.  Smagorinsky)  may 
be  applicable.  However,  such  standard  models  are  not  likely  to  be  able  to  produce  the  oscillatory 
behavior  in  the  dissipation  (figure  29),  as  our  experience  with  dynamic  Smagorinsky  in  this  context 
shows.  This  suggests  that  the  Fourier  cut-off  filter  would  be  a  poor  choice  for  wall  filtering,  if 
standard  subgrid  models  are  to  be  used. 

Finally,  we  note  that  the  tests  performed  here  were  particularly  arduous  for  the  filtered  bound¬ 
ary  formulation  because  discontinuities  in  derivatives  (as  in  the  velocity  at  the  wall)  are  poorly 
represented  by  Fourier  spectral  methods,  with  Gibbs  phenomena  as  the  result.  It  is  particularly 
remarkable,  then,  that  the  wall  stress  model  used  here  is  able  to  treat  and  largely  cancel  this  Gibbs 
phenomenon. 
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